src/gromacs/mdlib/forcerec.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "forcerec.h"
  40
  41 #include "config.h"
  42
  43 #include <cassert>
  44 #include <cmath>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include <algorithm>
  49 #include <memory>
  50
  51 #include "gromacs/commandline/filenm.h"
  52 #include "gromacs/domdec/domdec.h"
  53 #include "gromacs/domdec/domdec_struct.h"
  54 #include "gromacs/ewald/ewald.h"
  55 #include "gromacs/ewald/ewald_utils.h"
  56 #include "gromacs/fileio/filetypes.h"
  57 #include "gromacs/gmxlib/network.h"
  58 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  59 #include "gromacs/gpu_utils/gpu_utils.h"
  60 #include "gromacs/hardware/hw_info.h"
  61 #include "gromacs/listed_forces/gpubonded.h"
  62 #include "gromacs/listed_forces/manage_threading.h"
  63 #include "gromacs/listed_forces/pairs.h"
  64 #include "gromacs/math/functions.h"
  65 #include "gromacs/math/units.h"
  66 #include "gromacs/math/utilities.h"
  67 #include "gromacs/math/vec.h"
  68 #include "gromacs/mdlib/dispersioncorrection.h"
  69 #include "gromacs/mdlib/force.h"
  70 #include "gromacs/mdlib/forcerec_threading.h"
  71 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  72 #include "gromacs/mdlib/md_support.h"
  73 #include "gromacs/mdlib/ns.h"
  74 #include "gromacs/mdlib/qmmm.h"
  75 #include "gromacs/mdlib/rf_util.h"
  76 #include "gromacs/mdlib/wall.h"
  77 #include "gromacs/mdtypes/commrec.h"
  78 #include "gromacs/mdtypes/fcdata.h"
  79 #include "gromacs/mdtypes/group.h"
  80 #include "gromacs/mdtypes/iforceprovider.h"
  81 #include "gromacs/mdtypes/inputrec.h"
  82 #include "gromacs/mdtypes/md_enums.h"
  83 #include "gromacs/nbnxm/gpu_data_mgmt.h"
  84 #include "gromacs/nbnxm/nbnxm.h"
  85 #include "gromacs/nbnxm/nbnxm_geometry.h"
  86 #include "gromacs/pbcutil/ishift.h"
  87 #include "gromacs/pbcutil/pbc.h"
  88 #include "gromacs/tables/forcetable.h"
  89 #include "gromacs/topology/mtop_util.h"
  90 #include "gromacs/trajectory/trajectoryframe.h"
  91 #include "gromacs/utility/cstringutil.h"
  92 #include "gromacs/utility/exceptions.h"
  93 #include "gromacs/utility/fatalerror.h"
  94 #include "gromacs/utility/gmxassert.h"
  95 #include "gromacs/utility/logger.h"
  96 #include "gromacs/utility/physicalnodecommunicator.h"
  97 #include "gromacs/utility/pleasecite.h"
  98 #include "gromacs/utility/smalloc.h"
  99 #include "gromacs/utility/strconvert.h"
 100
 101 static real *mk_nbfp(const gmx_ffparams_t *idef, gmx_bool bBHAM)
 102 {
 103     real *nbfp;
 104     int   i, j, k, atnr;
 105
 106     atnr = idef->atnr;
 107     if (bBHAM)
 108     {
 109         snew(nbfp, 3*atnr*atnr);
 110         for (i = k = 0; (i < atnr); i++)
 111         {
 112             for (j = 0; (j < atnr); j++, k++)
 113             {
 114                 BHAMA(nbfp, atnr, i, j) = idef->iparams[k].bham.a;
 115                 BHAMB(nbfp, atnr, i, j) = idef->iparams[k].bham.b;
 116                 /* nbfp now includes the 6.0 derivative prefactor */
 117                 BHAMC(nbfp, atnr, i, j) = idef->iparams[k].bham.c*6.0;
 118             }
 119         }
 120     }
 121     else
 122     {
 123         snew(nbfp, 2*atnr*atnr);
 124         for (i = k = 0; (i < atnr); i++)
 125         {
 126             for (j = 0; (j < atnr); j++, k++)
 127             {
 128                 /* nbfp now includes the 6.0/12.0 derivative prefactors */
 129                 C6(nbfp, atnr, i, j)   = idef->iparams[k].lj.c6*6.0;
 130                 C12(nbfp, atnr, i, j)  = idef->iparams[k].lj.c12*12.0;
 131             }
 132         }
 133     }
 134
 135     return nbfp;
 136 }
 137
 138 static real *make_ljpme_c6grid(const gmx_ffparams_t *idef, t_forcerec *fr)
 139 {
 140     int        i, j, k, atnr;
 141     real       c6, c6i, c6j, c12i, c12j, epsi, epsj, sigmai, sigmaj;
 142     real      *grid;
 143
 144     /* For LJ-PME simulations, we correct the energies with the reciprocal space
 145      * inside of the cut-off. To do this the non-bonded kernels needs to have
 146      * access to the C6-values used on the reciprocal grid in pme.c
 147      */
 148
 149     atnr = idef->atnr;
 150     snew(grid, 2*atnr*atnr);
 151     for (i = k = 0; (i < atnr); i++)
 152     {
 153         for (j = 0; (j < atnr); j++, k++)
 154         {
 155             c6i  = idef->iparams[i*(atnr+1)].lj.c6;
 156             c12i = idef->iparams[i*(atnr+1)].lj.c12;
 157             c6j  = idef->iparams[j*(atnr+1)].lj.c6;
 158             c12j = idef->iparams[j*(atnr+1)].lj.c12;
 159             c6   = std::sqrt(c6i * c6j);
 160             if (fr->ljpme_combination_rule == eljpmeLB
 161                 && !gmx_numzero(c6) && !gmx_numzero(c12i) && !gmx_numzero(c12j))
 162             {
 163                 sigmai = gmx::sixthroot(c12i / c6i);
 164                 sigmaj = gmx::sixthroot(c12j / c6j);
 165                 epsi   = c6i * c6i / c12i;
 166                 epsj   = c6j * c6j / c12j;
 167                 c6     = std::sqrt(epsi * epsj) * gmx::power6(0.5*(sigmai+sigmaj));
 168             }
 169             /* Store the elements at the same relative positions as C6 in nbfp in order
 170              * to simplify access in the kernels
 171              */
 172             grid[2*(atnr*i+j)] = c6*6.0;
 173         }
 174     }
 175     return grid;
 176 }
 177
 178 /* This routine sets fr->solvent_opt to the most common solvent in the
 179  * system, e.g. esolSPC or esolTIP4P. It will also mark each charge group in
 180  * the fr->solvent_type array with the correct type (or esolNO).
 181  *
 182  * Charge groups that fulfill the conditions but are not identical to the
 183  * most common one will be marked as esolNO in the solvent_type array.
 184  *
 185  * TIP3p is identical to SPC for these purposes, so we call it
 186  * SPC in the arrays (Apologies to Bill Jorgensen ;-)
 187  *
 188  * NOTE: QM particle should not
 189  * become an optimized solvent. Not even if there is only one charge
 190  * group in the Qm
 191  */
 192
 193 typedef struct
 194 {
 195     int    model;
 196     int    count;
 197     int    vdwtype[4];
 198     real   charge[4];
 199 } solvent_parameters_t;
 200
 201 static void
 202 check_solvent_cg(const gmx_moltype_t    *molt,
 203                  int                     cg0,
 204                  int                     nmol,
 205                  const unsigned char    *qm_grpnr,
 206                  const t_grps           *qm_grps,
 207                  t_forcerec   *          fr,
 208                  int                    *n_solvent_parameters,
 209                  solvent_parameters_t  **solvent_parameters_p,
 210                  int                     cginfo,
 211                  int                    *cg_sp)
 212 {
 213     t_atom               *atom;
 214     int                   j, k;
 215     int                   j0, j1, nj;
 216     gmx_bool              perturbed;
 217     gmx_bool              has_vdw[4];
 218     gmx_bool              match;
 219     real                  tmp_charge[4]  = { 0.0 }; /* init to zero to make gcc 7 happy */
 220     int                   tmp_vdwtype[4] = { 0 };   /* init to zero to make gcc 7 happy */
 221     int                   tjA;
 222     gmx_bool              qm;
 223     solvent_parameters_t *solvent_parameters;
 224
 225     /* We use a list with parameters for each solvent type.
 226      * Every time we discover a new molecule that fulfills the basic
 227      * conditions for a solvent we compare with the previous entries
 228      * in these lists. If the parameters are the same we just increment
 229      * the counter for that type, and otherwise we create a new type
 230      * based on the current molecule.
 231      *
 232      * Once we've finished going through all molecules we check which
 233      * solvent is most common, and mark all those molecules while we
 234      * clear the flag on all others.
 235      */
 236
 237     solvent_parameters = *solvent_parameters_p;
 238
 239     /* Mark the cg first as non optimized */
 240     *cg_sp = -1;
 241
 242     /* Check if this cg has no exclusions with atoms in other charge groups
 243      * and all atoms inside the charge group excluded.
 244      * We only have 3 or 4 atom solvent loops.
 245      */
 246     if (GET_CGINFO_EXCL_INTER(cginfo) ||
 247         !GET_CGINFO_EXCL_INTRA(cginfo))
 248     {
 249         return;
 250     }
 251
 252     /* Get the indices of the first atom in this charge group */
 253     j0     = molt->cgs.index[cg0];
 254     j1     = molt->cgs.index[cg0+1];
 255
 256     /* Number of atoms in our molecule */
 257     nj     = j1 - j0;
 258
 259     if (debug)
 260     {
 261         fprintf(debug,
 262                 "Moltype '%s': there are %d atoms in this charge group\n",
 263                 *molt->name, nj);
 264     }
 265
 266     /* Check if it could be an SPC (3 atoms) or TIP4p (4) water,
 267      * otherwise skip it.
 268      */
 269     if (nj < 3 || nj > 4)
 270     {
 271         return;
 272     }
 273
 274     /* Check if we are doing QM on this group */
 275     qm = FALSE;
 276     if (qm_grpnr != nullptr)
 277     {
 278         for (j = j0; j < j1 && !qm; j++)
 279         {
 280             qm = (qm_grpnr[j] < qm_grps->nr - 1);
 281         }
 282     }
 283     /* Cannot use solvent optimization with QM */
 284     if (qm)
 285     {
 286         return;
 287     }
 288
 289     atom = molt->atoms.atom;
 290
 291     /* Still looks like a solvent, time to check parameters */
 292
 293     /* If it is perturbed (free energy) we can't use the solvent loops,
 294      * so then we just skip to the next molecule.
 295      */
 296     perturbed = FALSE;
 297
 298     for (j = j0; j < j1 && !perturbed; j++)
 299     {
 300         perturbed = PERTURBED(atom[j]);
 301     }
 302
 303     if (perturbed)
 304     {
 305         return;
 306     }
 307
 308     /* Now it's only a question if the VdW and charge parameters
 309      * are OK. Before doing the check we compare and see if they are
 310      * identical to a possible previous solvent type.
 311      * First we assign the current types and charges.
 312      */
 313     for (j = 0; j < nj; j++)
 314     {
 315         tmp_vdwtype[j] = atom[j0+j].type;
 316         tmp_charge[j]  = atom[j0+j].q;
 317     }
 318
 319     /* Does it match any previous solvent type? */
 320     for (k = 0; k < *n_solvent_parameters; k++)
 321     {
 322         match = TRUE;
 323
 324
 325         /* We can only match SPC with 3 atoms and TIP4p with 4 atoms */
 326         if ( (solvent_parameters[k].model == esolSPC   && nj != 3)  ||
 327              (solvent_parameters[k].model == esolTIP4P && nj != 4) )
 328         {
 329             match = FALSE;
 330         }
 331
 332         /* Check that types & charges match for all atoms in molecule */
 333         for (j = 0; j < nj && match; j++)
 334         {
 335             if (tmp_vdwtype[j] != solvent_parameters[k].vdwtype[j])
 336             {
 337                 match = FALSE;
 338             }
 339             if (tmp_charge[j] != solvent_parameters[k].charge[j])
 340             {
 341                 match = FALSE;
 342             }
 343         }
 344         if (match)
 345         {
 346             /* Congratulations! We have a matched solvent.
 347              * Flag it with this type for later processing.
 348              */
 349             *cg_sp = k;
 350             solvent_parameters[k].count += nmol;
 351
 352             /* We are done with this charge group */
 353             return;
 354         }
 355     }
 356
 357     /* If we get here, we have a tentative new solvent type.
 358      * Before we add it we must check that it fulfills the requirements
 359      * of the solvent optimized loops. First determine which atoms have
 360      * VdW interactions.
 361      */
 362     for (j = 0; j < nj; j++)
 363     {
 364         has_vdw[j] = FALSE;
 365         tjA        = tmp_vdwtype[j];
 366
 367         /* Go through all other tpes and see if any have non-zero
 368          * VdW parameters when combined with this one.
 369          */
 370         for (k = 0; k < fr->ntype && (!has_vdw[j]); k++)
 371         {
 372             /* We already checked that the atoms weren't perturbed,
 373              * so we only need to check state A now.
 374              */
 375             if (fr->bBHAM)
 376             {
 377                 has_vdw[j] = (has_vdw[j] ||
 378                               (BHAMA(fr->nbfp, fr->ntype, tjA, k) != 0.0) ||
 379                               (BHAMB(fr->nbfp, fr->ntype, tjA, k) != 0.0) ||
 380                               (BHAMC(fr->nbfp, fr->ntype, tjA, k) != 0.0));
 381             }
 382             else
 383             {
 384                 /* Standard LJ */
 385                 has_vdw[j] = (has_vdw[j] ||
 386                               (C6(fr->nbfp, fr->ntype, tjA, k)  != 0.0) ||
 387                               (C12(fr->nbfp, fr->ntype, tjA, k) != 0.0));
 388             }
 389         }
 390     }
 391
 392     /* Now we know all we need to make the final check and assignment. */
 393     if (nj == 3)
 394     {
 395         /* So, is it an SPC?
 396          * For this we require thatn all atoms have charge,
 397          * the charges on atom 2 & 3 should be the same, and only
 398          * atom 1 might have VdW.
 399          */
 400         if (!has_vdw[1] &&
 401             !has_vdw[2] &&
 402             tmp_charge[0]  != 0 &&
 403             tmp_charge[1]  != 0 &&
 404             tmp_charge[2]  == tmp_charge[1])
 405         {
 406             srenew(solvent_parameters, *n_solvent_parameters+1);
 407             solvent_parameters[*n_solvent_parameters].model = esolSPC;
 408             solvent_parameters[*n_solvent_parameters].count = nmol;
 409             for (k = 0; k < 3; k++)
 410             {
 411                 solvent_parameters[*n_solvent_parameters].vdwtype[k] = tmp_vdwtype[k];
 412                 solvent_parameters[*n_solvent_parameters].charge[k]  = tmp_charge[k];
 413             }
 414
 415             *cg_sp = *n_solvent_parameters;
 416             (*n_solvent_parameters)++;
 417         }
 418     }
 419     else if (nj == 4)
 420     {
 421         /* Or could it be a TIP4P?
 422          * For this we require thatn atoms 2,3,4 have charge, but not atom 1.
 423          * Only atom 1 mght have VdW.
 424          */
 425         if (!has_vdw[1] &&
 426             !has_vdw[2] &&
 427             !has_vdw[3] &&
 428             tmp_charge[0]  == 0 &&
 429             tmp_charge[1]  != 0 &&
 430             tmp_charge[2]  == tmp_charge[1] &&
 431             tmp_charge[3]  != 0)
 432         {
 433             srenew(solvent_parameters, *n_solvent_parameters+1);
 434             solvent_parameters[*n_solvent_parameters].model = esolTIP4P;
 435             solvent_parameters[*n_solvent_parameters].count = nmol;
 436             for (k = 0; k < 4; k++)
 437             {
 438                 solvent_parameters[*n_solvent_parameters].vdwtype[k] = tmp_vdwtype[k];
 439                 solvent_parameters[*n_solvent_parameters].charge[k]  = tmp_charge[k];
 440             }
 441
 442             *cg_sp = *n_solvent_parameters;
 443             (*n_solvent_parameters)++;
 444         }
 445     }
 446
 447     *solvent_parameters_p = solvent_parameters;
 448 }
 449
 450 static void
 451 check_solvent(FILE  *                fp,
 452               const gmx_mtop_t  *    mtop,
 453               t_forcerec  *          fr,
 454               cginfo_mb_t           *cginfo_mb)
 455 {
 456     const t_block     *   cgs;
 457     const gmx_moltype_t  *molt;
 458     int                   mol, cg_mol, at_offset, am, cgm, i, nmol_ch, nmol;
 459     int                   n_solvent_parameters;
 460     solvent_parameters_t *solvent_parameters;
 461     int                 **cg_sp;
 462     int                   bestsp, bestsol;
 463
 464     if (debug)
 465     {
 466         fprintf(debug, "Going to determine what solvent types we have.\n");
 467     }
 468
 469     n_solvent_parameters = 0;
 470     solvent_parameters   = nullptr;
 471     /* Allocate temporary array for solvent type */
 472     snew(cg_sp, mtop->molblock.size());
 473
 474     at_offset = 0;
 475     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 476     {
 477         molt = &mtop->moltype[mtop->molblock[mb].type];
 478         cgs  = &molt->cgs;
 479         /* Here we have to loop over all individual molecules
 480          * because we need to check for QMMM particles.
 481          */
 482         snew(cg_sp[mb], cginfo_mb[mb].cg_mod);
 483         nmol_ch = cginfo_mb[mb].cg_mod/cgs->nr;
 484         nmol    = mtop->molblock[mb].nmol/nmol_ch;
 485         for (mol = 0; mol < nmol_ch; mol++)
 486         {
 487             cgm = mol*cgs->nr;
 488             am  = mol*cgs->index[cgs->nr];
 489             for (cg_mol = 0; cg_mol < cgs->nr; cg_mol++)
 490             {
 491                 check_solvent_cg(molt, cg_mol, nmol,
 492                                  mtop->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].empty() ?
 493                                  nullptr : mtop->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].data()+at_offset+am,
 494                                  &mtop->groups.groups[SimulationAtomGroupType::QuantumMechanics],
 495                                  fr,
 496                                  &n_solvent_parameters, &solvent_parameters,
 497                                  cginfo_mb[mb].cginfo[cgm+cg_mol],
 498                                  &cg_sp[mb][cgm+cg_mol]);
 499             }
 500         }
 501         at_offset += cgs->index[cgs->nr];
 502     }
 503
 504     /* Puh! We finished going through all charge groups.
 505      * Now find the most common solvent model.
 506      */
 507
 508     /* Most common solvent this far */
 509     bestsp = -2;
 510     for (i = 0; i < n_solvent_parameters; i++)
 511     {
 512         if (bestsp == -2 ||
 513             solvent_parameters[i].count > solvent_parameters[bestsp].count)
 514         {
 515             bestsp = i;
 516         }
 517     }
 518
 519     if (bestsp >= 0)
 520     {
 521         bestsol = solvent_parameters[bestsp].model;
 522     }
 523     else
 524     {
 525         bestsol = esolNO;
 526     }
 527
 528     fr->nWatMol = 0;
 529     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 530     {
 531         cgs  = &mtop->moltype[mtop->molblock[mb].type].cgs;
 532         nmol = (mtop->molblock[mb].nmol*cgs->nr)/cginfo_mb[mb].cg_mod;
 533         for (i = 0; i < cginfo_mb[mb].cg_mod; i++)
 534         {
 535             if (cg_sp[mb][i] == bestsp)
 536             {
 537                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[i], bestsol);
 538                 fr->nWatMol += nmol;
 539             }
 540             else
 541             {
 542                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[i], esolNO);
 543             }
 544         }
 545         sfree(cg_sp[mb]);
 546     }
 547     sfree(cg_sp);
 548
 549     if (bestsol != esolNO && fp != nullptr)
 550     {
 551         fprintf(fp, "\nEnabling %s-like water optimization for %d molecules.\n\n",
 552                 esol_names[bestsol],
 553                 solvent_parameters[bestsp].count);
 554     }
 555
 556     sfree(solvent_parameters);
 557     fr->solvent_opt = bestsol;
 558 }
 559
 560 enum {
 561     acNONE = 0, acCONSTRAINT, acSETTLE
 562 };
 563
 564 static cginfo_mb_t *init_cginfo_mb(FILE *fplog, const gmx_mtop_t *mtop,
 565                                    t_forcerec *fr, gmx_bool bNoSolvOpt,
 566                                    gmx_bool *bFEP_NonBonded,
 567                                    gmx_bool *bExcl_IntraCGAll_InterCGNone)
 568 {
 569     const t_block        *cgs;
 570     const t_blocka       *excl;
 571     const gmx_moltype_t  *molt;
 572     const gmx_molblock_t *molb;
 573     cginfo_mb_t          *cginfo_mb;
 574     gmx_bool             *type_VDW;
 575     int                  *cginfo;
 576     int                   cg_offset, a_offset;
 577     int                   m, cg, a0, a1, gid, ai, j, aj, excl_nalloc;
 578     int                  *a_con;
 579     int                   ftype;
 580     int                   ia;
 581     gmx_bool              bId, *bExcl, bExclIntraAll, bExclInter, bHaveVDW, bHaveQ, bHavePerturbedAtoms;
 582
 583     snew(cginfo_mb, mtop->molblock.size());
 584
 585     snew(type_VDW, fr->ntype);
 586     for (ai = 0; ai < fr->ntype; ai++)
 587     {
 588         type_VDW[ai] = FALSE;
 589         for (j = 0; j < fr->ntype; j++)
 590         {
 591             type_VDW[ai] = type_VDW[ai] ||
 592                 fr->bBHAM ||
 593                 C6(fr->nbfp, fr->ntype, ai, j) != 0 ||
 594                 C12(fr->nbfp, fr->ntype, ai, j) != 0;
 595         }
 596     }
 597
 598     *bFEP_NonBonded               = FALSE;
 599     *bExcl_IntraCGAll_InterCGNone = TRUE;
 600
 601     excl_nalloc = 10;
 602     snew(bExcl, excl_nalloc);
 603     cg_offset = 0;
 604     a_offset  = 0;
 605     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 606     {
 607         molb = &mtop->molblock[mb];
 608         molt = &mtop->moltype[molb->type];
 609         cgs  = &molt->cgs;
 610         excl = &molt->excls;
 611
 612         /* Check if the cginfo is identical for all molecules in this block.
 613          * If so, we only need an array of the size of one molecule.
 614          * Otherwise we make an array of #mol times #cgs per molecule.
 615          */
 616         bId = TRUE;
 617         for (m = 0; m < molb->nmol; m++)
 618         {
 619             int am = m*cgs->index[cgs->nr];
 620             for (cg = 0; cg < cgs->nr; cg++)
 621             {
 622                 a0 = cgs->index[cg];
 623                 a1 = cgs->index[cg+1];
 624                 if (getGroupType(mtop->groups, SimulationAtomGroupType::QuantumMechanics, a_offset+am+a0) !=
 625                     getGroupType(mtop->groups, SimulationAtomGroupType::QuantumMechanics, a_offset   +a0))
 626                 {
 627                     bId = FALSE;
 628                 }
 629                 if (!mtop->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].empty())
 630                 {
 631                     for (ai = a0; ai < a1; ai++)
 632                     {
 633                         if (mtop->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics][a_offset+am+ai] !=
 634                             mtop->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics][a_offset   +ai])
 635                         {
 636                             bId = FALSE;
 637                         }
 638                     }
 639                 }
 640             }
 641         }
 642
 643         cginfo_mb[mb].cg_start = cg_offset;
 644         cginfo_mb[mb].cg_end   = cg_offset + molb->nmol*cgs->nr;
 645         cginfo_mb[mb].cg_mod   = (bId ? 1 : molb->nmol)*cgs->nr;
 646         snew(cginfo_mb[mb].cginfo, cginfo_mb[mb].cg_mod);
 647         cginfo = cginfo_mb[mb].cginfo;
 648
 649         /* Set constraints flags for constrained atoms */
 650         snew(a_con, molt->atoms.nr);
 651         for (ftype = 0; ftype < F_NRE; ftype++)
 652         {
 653             if (interaction_function[ftype].flags & IF_CONSTRAINT)
 654             {
 655                 int nral;
 656
 657                 nral = NRAL(ftype);
 658                 for (ia = 0; ia < molt->ilist[ftype].size(); ia += 1+nral)
 659                 {
 660                     int a;
 661
 662                     for (a = 0; a < nral; a++)
 663                     {
 664                         a_con[molt->ilist[ftype].iatoms[ia+1+a]] =
 665                             (ftype == F_SETTLE ? acSETTLE : acCONSTRAINT);
 666                     }
 667                 }
 668             }
 669         }
 670
 671         for (m = 0; m < (bId ? 1 : molb->nmol); m++)
 672         {
 673             int cgm = m*cgs->nr;
 674             int am  = m*cgs->index[cgs->nr];
 675             for (cg = 0; cg < cgs->nr; cg++)
 676             {
 677                 a0 = cgs->index[cg];
 678                 a1 = cgs->index[cg+1];
 679
 680                 /* Store the energy group in cginfo */
 681                 gid = getGroupType(mtop->groups, SimulationAtomGroupType::EnergyOutput, a_offset+am+a0);
 682                 SET_CGINFO_GID(cginfo[cgm+cg], gid);
 683
 684                 /* Check the intra/inter charge group exclusions */
 685                 if (a1-a0 > excl_nalloc)
 686                 {
 687                     excl_nalloc = a1 - a0;
 688                     srenew(bExcl, excl_nalloc);
 689                 }
 690                 /* bExclIntraAll: all intra cg interactions excluded
 691                  * bExclInter:    any inter cg interactions excluded
 692                  */
 693                 bExclIntraAll       = TRUE;
 694                 bExclInter          = FALSE;
 695                 bHaveVDW            = FALSE;
 696                 bHaveQ              = FALSE;
 697                 bHavePerturbedAtoms = FALSE;
 698                 for (ai = a0; ai < a1; ai++)
 699                 {
 700                     /* Check VDW and electrostatic interactions */
 701                     bHaveVDW = bHaveVDW || (type_VDW[molt->atoms.atom[ai].type] ||
 702                                             type_VDW[molt->atoms.atom[ai].typeB]);
 703                     bHaveQ  = bHaveQ    || (molt->atoms.atom[ai].q != 0 ||
 704                                             molt->atoms.atom[ai].qB != 0);
 705
 706                     bHavePerturbedAtoms = bHavePerturbedAtoms || (PERTURBED(molt->atoms.atom[ai]) != 0);
 707
 708                     /* Clear the exclusion list for atom ai */
 709                     for (aj = a0; aj < a1; aj++)
 710                     {
 711                         bExcl[aj-a0] = FALSE;
 712                     }
 713                     /* Loop over all the exclusions of atom ai */
 714                     for (j = excl->index[ai]; j < excl->index[ai+1]; j++)
 715                     {
 716                         aj = excl->a[j];
 717                         if (aj < a0 || aj >= a1)
 718                         {
 719                             bExclInter = TRUE;
 720                         }
 721                         else
 722                         {
 723                             bExcl[aj-a0] = TRUE;
 724                         }
 725                     }
 726                     /* Check if ai excludes a0 to a1 */
 727                     for (aj = a0; aj < a1; aj++)
 728                     {
 729                         if (!bExcl[aj-a0])
 730                         {
 731                             bExclIntraAll = FALSE;
 732                         }
 733                     }
 734
 735                     switch (a_con[ai])
 736                     {
 737                         case acCONSTRAINT:
 738                             SET_CGINFO_CONSTR(cginfo[cgm+cg]);
 739                             break;
 740                         case acSETTLE:
 741                             SET_CGINFO_SETTLE(cginfo[cgm+cg]);
 742                             break;
 743                         default:
 744                             break;
 745                     }
 746                 }
 747                 if (bExclIntraAll)
 748                 {
 749                     SET_CGINFO_EXCL_INTRA(cginfo[cgm+cg]);
 750                 }
 751                 if (bExclInter)
 752                 {
 753                     SET_CGINFO_EXCL_INTER(cginfo[cgm+cg]);
 754                 }
 755                 if (a1 - a0 > MAX_CHARGEGROUP_SIZE)
 756                 {
 757                     /* The size in cginfo is currently only read with DD */
 758                     gmx_fatal(FARGS, "A charge group has size %d which is larger than the limit of %d atoms", a1-a0, MAX_CHARGEGROUP_SIZE);
 759                 }
 760                 if (bHaveVDW)
 761                 {
 762                     SET_CGINFO_HAS_VDW(cginfo[cgm+cg]);
 763                 }
 764                 if (bHaveQ)
 765                 {
 766                     SET_CGINFO_HAS_Q(cginfo[cgm+cg]);
 767                 }
 768                 if (bHavePerturbedAtoms && fr->efep != efepNO)
 769                 {
 770                     SET_CGINFO_FEP(cginfo[cgm+cg]);
 771                     *bFEP_NonBonded = TRUE;
 772                 }
 773                 /* Store the charge group size */
 774                 SET_CGINFO_NATOMS(cginfo[cgm+cg], a1-a0);
 775
 776                 if (!bExclIntraAll || bExclInter)
 777                 {
 778                     *bExcl_IntraCGAll_InterCGNone = FALSE;
 779                 }
 780             }
 781         }
 782
 783         sfree(a_con);
 784
 785         cg_offset += molb->nmol*cgs->nr;
 786         a_offset  += molb->nmol*cgs->index[cgs->nr];
 787     }
 788     sfree(type_VDW);
 789     sfree(bExcl);
 790
 791     /* the solvent optimizer is called after the QM is initialized,
 792      * because we don't want to have the QM subsystemto become an
 793      * optimized solvent
 794      */
 795
 796     check_solvent(fplog, mtop, fr, cginfo_mb);
 797
 798     if (getenv("GMX_NO_SOLV_OPT"))
 799     {
 800         if (fplog)
 801         {
 802             fprintf(fplog, "Found environment variable GMX_NO_SOLV_OPT.\n"
 803                     "Disabling all solvent optimization\n");
 804         }
 805         fr->solvent_opt = esolNO;
 806     }
 807     if (bNoSolvOpt)
 808     {
 809         fr->solvent_opt = esolNO;
 810     }
 811     if (!fr->solvent_opt)
 812     {
 813         for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 814         {
 815             for (cg = 0; cg < cginfo_mb[mb].cg_mod; cg++)
 816             {
 817                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[cg], esolNO);
 818             }
 819         }
 820     }
 821
 822     return cginfo_mb;
 823 }
 824
 825 static int *cginfo_expand(int nmb, cginfo_mb_t *cgi_mb)
 826 {
 827     int  ncg, mb, cg;
 828     int *cginfo;
 829
 830     ncg = cgi_mb[nmb-1].cg_end;
 831     snew(cginfo, ncg);
 832     mb = 0;
 833     for (cg = 0; cg < ncg; cg++)
 834     {
 835         while (cg >= cgi_mb[mb].cg_end)
 836         {
 837             mb++;
 838         }
 839         cginfo[cg] =
 840             cgi_mb[mb].cginfo[(cg - cgi_mb[mb].cg_start) % cgi_mb[mb].cg_mod];
 841     }
 842
 843     return cginfo;
 844 }
 845
 846 static void done_cginfo_mb(cginfo_mb_t *cginfo_mb, int numMolBlocks)
 847 {
 848     if (cginfo_mb == nullptr)
 849     {
 850         return;
 851     }
 852     for (int mb = 0; mb < numMolBlocks; ++mb)
 853     {
 854         sfree(cginfo_mb[mb].cginfo);
 855     }
 856     sfree(cginfo_mb);
 857 }
 858
 859 /* Sets the sum of charges (squared) and C6 in the system in fr.
 860  * Returns whether the system has a net charge.
 861  */
 862 static bool set_chargesum(FILE *log, t_forcerec *fr, const gmx_mtop_t *mtop)
 863 {
 864     /*This now calculates sum for q and c6*/
 865     double         qsum, q2sum, q, c6sum, c6;
 866
 867     qsum   = 0;
 868     q2sum  = 0;
 869     c6sum  = 0;
 870     for (const gmx_molblock_t &molb : mtop->molblock)
 871     {
 872         int            nmol  = molb.nmol;
 873         const t_atoms *atoms = &mtop->moltype[molb.type].atoms;
 874         for (int i = 0; i < atoms->nr; i++)
 875         {
 876             q       = atoms->atom[i].q;
 877             qsum   += nmol*q;
 878             q2sum  += nmol*q*q;
 879             c6      = mtop->ffparams.iparams[atoms->atom[i].type*(mtop->ffparams.atnr+1)].lj.c6;
 880             c6sum  += nmol*c6;
 881         }
 882     }
 883     fr->qsum[0]   = qsum;
 884     fr->q2sum[0]  = q2sum;
 885     fr->c6sum[0]  = c6sum;
 886
 887     if (fr->efep != efepNO)
 888     {
 889         qsum   = 0;
 890         q2sum  = 0;
 891         c6sum  = 0;
 892         for (const gmx_molblock_t &molb : mtop->molblock)
 893         {
 894             int            nmol  = molb.nmol;
 895             const t_atoms *atoms = &mtop->moltype[molb.type].atoms;
 896             for (int i = 0; i < atoms->nr; i++)
 897             {
 898                 q       = atoms->atom[i].qB;
 899                 qsum   += nmol*q;
 900                 q2sum  += nmol*q*q;
 901                 c6      = mtop->ffparams.iparams[atoms->atom[i].typeB*(mtop->ffparams.atnr+1)].lj.c6;
 902                 c6sum  += nmol*c6;
 903             }
 904             fr->qsum[1]   = qsum;
 905             fr->q2sum[1]  = q2sum;
 906             fr->c6sum[1]  = c6sum;
 907         }
 908     }
 909     else
 910     {
 911         fr->qsum[1]   = fr->qsum[0];
 912         fr->q2sum[1]  = fr->q2sum[0];
 913         fr->c6sum[1]  = fr->c6sum[0];
 914     }
 915     if (log)
 916     {
 917         if (fr->efep == efepNO)
 918         {
 919             fprintf(log, "System total charge: %.3f\n", fr->qsum[0]);
 920         }
 921         else
 922         {
 923             fprintf(log, "System total charge, top. A: %.3f top. B: %.3f\n",
 924                     fr->qsum[0], fr->qsum[1]);
 925         }
 926     }
 927
 928     /* A cut-off of 1e-4 is used to catch rounding errors due to ascii input */
 929     return (std::abs(fr->qsum[0]) > 1e-4 ||
 930             std::abs(fr->qsum[1]) > 1e-4);
 931 }
 932
 933 void update_forcerec(t_forcerec *fr, matrix box)
 934 {
 935     if (fr->ic->eeltype == eelGRF)
 936     {
 937         calc_rffac(nullptr, fr->ic->eeltype, fr->ic->epsilon_r, fr->ic->epsilon_rf,
 938                    fr->ic->rcoulomb, fr->temp, fr->zsquare, box,
 939                    &fr->ic->k_rf, &fr->ic->c_rf);
 940     }
 941 }
 942
 943 static real calcBuckinghamBMax(FILE *fplog, const gmx_mtop_t *mtop)
 944 {
 945     const t_atoms *at1, *at2;
 946     int            i, j, tpi, tpj, ntypes;
 947     real           b, bmin;
 948
 949     if (fplog)
 950     {
 951         fprintf(fplog, "Determining largest Buckingham b parameter for table\n");
 952     }
 953     ntypes = mtop->ffparams.atnr;
 954
 955     bmin            = -1;
 956     real bham_b_max = 0;
 957     for (size_t mt1 = 0; mt1 < mtop->moltype.size(); mt1++)
 958     {
 959         at1 = &mtop->moltype[mt1].atoms;
 960         for (i = 0; (i < at1->nr); i++)
 961         {
 962             tpi = at1->atom[i].type;
 963             if (tpi >= ntypes)
 964             {
 965                 gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", i, tpi, ntypes);
 966             }
 967
 968             for (size_t mt2 = mt1; mt2 < mtop->moltype.size(); mt2++)
 969             {
 970                 at2 = &mtop->moltype[mt2].atoms;
 971                 for (j = 0; (j < at2->nr); j++)
 972                 {
 973                     tpj = at2->atom[j].type;
 974                     if (tpj >= ntypes)
 975                     {
 976                         gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", j, tpj, ntypes);
 977                     }
 978                     b = mtop->ffparams.iparams[tpi*ntypes + tpj].bham.b;
 979                     if (b > bham_b_max)
 980                     {
 981                         bham_b_max = b;
 982                     }
 983                     if ((b < bmin) || (bmin == -1))
 984                     {
 985                         bmin = b;
 986                     }
 987                 }
 988             }
 989         }
 990     }
 991     if (fplog)
 992     {
 993         fprintf(fplog, "Buckingham b parameters, min: %g, max: %g\n",
 994                 bmin, bham_b_max);
 995     }
 996
 997     return bham_b_max;
 998 }
 999
1000 static void make_nbf_tables(FILE *fp,
1001                             const interaction_const_t *ic, real rtab,
1002                             const char *tabfn, char *eg1, char *eg2,
1003                             t_nblists *nbl)
1004 {
1005     char buf[STRLEN];
1006     int  i, j;
1007
1008     if (tabfn == nullptr)
1009     {
1010         if (debug)
1011         {
1012             fprintf(debug, "No table file name passed, can not read table, can not do non-bonded interactions\n");
1013         }
1014         return;
1015     }
1016
1017     sprintf(buf, "%s", tabfn);
1018     if (eg1 && eg2)
1019     {
1020         /* Append the two energy group names */
1021         sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1, "_%s_%s.%s",
1022                 eg1, eg2, ftp2ext(efXVG));
1023     }
1024     nbl->table_elec_vdw = make_tables(fp, ic, buf, rtab, 0);
1025     /* Copy the contents of the table to separate coulomb and LJ tables too,
1026      * to improve cache performance.
1027      */
1028     /* For performance reasons we want
1029      * the table data to be aligned to 16-byte. The pointers could be freed
1030      * but currently aren't.
1031      */
1032     snew(nbl->table_elec, 1);
1033     nbl->table_elec->interaction   = GMX_TABLE_INTERACTION_ELEC;
1034     nbl->table_elec->format        = nbl->table_elec_vdw->format;
1035     nbl->table_elec->r             = nbl->table_elec_vdw->r;
1036     nbl->table_elec->n             = nbl->table_elec_vdw->n;
1037     nbl->table_elec->scale         = nbl->table_elec_vdw->scale;
1038     nbl->table_elec->formatsize    = nbl->table_elec_vdw->formatsize;
1039     nbl->table_elec->ninteractions = 1;
1040     nbl->table_elec->stride        = nbl->table_elec->formatsize * nbl->table_elec->ninteractions;
1041     snew_aligned(nbl->table_elec->data, nbl->table_elec->stride*(nbl->table_elec->n+1), 32);
1042
1043     snew(nbl->table_vdw, 1);
1044     nbl->table_vdw->interaction   = GMX_TABLE_INTERACTION_VDWREP_VDWDISP;
1045     nbl->table_vdw->format        = nbl->table_elec_vdw->format;
1046     nbl->table_vdw->r             = nbl->table_elec_vdw->r;
1047     nbl->table_vdw->n             = nbl->table_elec_vdw->n;
1048     nbl->table_vdw->scale         = nbl->table_elec_vdw->scale;
1049     nbl->table_vdw->formatsize    = nbl->table_elec_vdw->formatsize;
1050     nbl->table_vdw->ninteractions = 2;
1051     nbl->table_vdw->stride        = nbl->table_vdw->formatsize * nbl->table_vdw->ninteractions;
1052     snew_aligned(nbl->table_vdw->data, nbl->table_vdw->stride*(nbl->table_vdw->n+1), 32);
1053
1054     /* NOTE: Using a single i-loop here leads to mix-up of data in table_vdw
1055      *       with (at least) gcc 6.2, 6.3 and 6.4 when compiled with -O3 and AVX
1056      */
1057     for (i = 0; i <= nbl->table_elec_vdw->n; i++)
1058     {
1059         for (j = 0; j < 4; j++)
1060         {
1061             nbl->table_elec->data[4*i+j] = nbl->table_elec_vdw->data[12*i+j];
1062         }
1063     }
1064     for (i = 0; i <= nbl->table_elec_vdw->n; i++)
1065     {
1066         for (j = 0; j < 8; j++)
1067         {
1068             nbl->table_vdw->data[8*i+j] = nbl->table_elec_vdw->data[12*i+4+j];
1069         }
1070     }
1071 }
1072
1073 /*!\brief If there's bonded interactions of type \c ftype1 or \c
1074  * ftype2 present in the topology, build an array of the number of
1075  * interactions present for each bonded interaction index found in the
1076  * topology.
1077  *
1078  * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
1079  * valid type with that parameter.
1080  *
1081  * \c count will be reallocated as necessary to fit the largest bonded
1082  * interaction index found, and its current size will be returned in
1083  * \c ncount. It will contain zero for every bonded interaction index
1084  * for which no interactions are present in the topology.
1085  */
1086 static void count_tables(int ftype1, int ftype2, const gmx_mtop_t *mtop,
1087                          int *ncount, int **count)
1088 {
1089     int                  ftype, i, j, tabnr;
1090
1091     // Loop over all moleculetypes
1092     for (const gmx_moltype_t &molt : mtop->moltype)
1093     {
1094         // Loop over all interaction types
1095         for (ftype = 0; ftype < F_NRE; ftype++)
1096         {
1097             // If the current interaction type is one of the types whose tables we're trying to count...
1098             if (ftype == ftype1 || ftype == ftype2)
1099             {
1100                 const InteractionList &il     = molt.ilist[ftype];
1101                 const int              stride = 1 + NRAL(ftype);
1102                 // ... and there are actually some interactions for this type
1103                 for (i = 0; i < il.size(); i += stride)
1104                 {
1105                     // Find out which table index the user wanted
1106                     tabnr = mtop->ffparams.iparams[il.iatoms[i]].tab.table;
1107                     if (tabnr < 0)
1108                     {
1109                         gmx_fatal(FARGS, "A bonded table number is smaller than 0: %d\n", tabnr);
1110                     }
1111                     // Make room for this index in the data structure
1112                     if (tabnr >= *ncount)
1113                     {
1114                         srenew(*count, tabnr+1);
1115                         for (j = *ncount; j < tabnr+1; j++)
1116                         {
1117                             (*count)[j] = 0;
1118                         }
1119                         *ncount = tabnr+1;
1120                     }
1121                     // Record that this table index is used and must have a valid file
1122                     (*count)[tabnr]++;
1123                 }
1124             }
1125         }
1126     }
1127 }
1128
1129 /*!\brief If there's bonded interactions of flavour \c tabext and type
1130  * \c ftype1 or \c ftype2 present in the topology, seek them in the
1131  * list of filenames passed to mdrun, and make bonded tables from
1132  * those files.
1133  *
1134  * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
1135  * valid type with that parameter.
1136  *
1137  * A fatal error occurs if no matching filename is found.
1138  */
1139 static bondedtable_t *make_bonded_tables(FILE *fplog,
1140                                          int ftype1, int ftype2,
1141                                          const gmx_mtop_t *mtop,
1142                                          gmx::ArrayRef<const std::string> tabbfnm,
1143                                          const char *tabext)
1144 {
1145     int            ncount, *count;
1146     bondedtable_t *tab;
1147
1148     tab = nullptr;
1149
1150     ncount = 0;
1151     count  = nullptr;
1152     count_tables(ftype1, ftype2, mtop, &ncount, &count);
1153
1154     // Are there any relevant tabulated bond interactions?
1155     if (ncount > 0)
1156     {
1157         snew(tab, ncount);
1158         for (int i = 0; i < ncount; i++)
1159         {
1160             // Do any interactions exist that requires this table?
1161             if (count[i] > 0)
1162             {
1163                 // This pattern enforces the current requirement that
1164                 // table filenames end in a characteristic sequence
1165                 // before the file type extension, and avoids table 13
1166                 // being recognized and used for table 1.
1167                 std::string patternToFind = gmx::formatString("_%s%d.%s", tabext, i, ftp2ext(efXVG));
1168                 bool        madeTable     = false;
1169                 for (gmx::index j = 0; j < tabbfnm.ssize() && !madeTable; ++j)
1170                 {
1171                     if (gmx::endsWith(tabbfnm[j], patternToFind))
1172                     {
1173                         // Finally read the table from the file found
1174                         tab[i]    = make_bonded_table(fplog, tabbfnm[j].c_str(), NRAL(ftype1)-2);
1175                         madeTable = true;
1176                     }
1177                 }
1178                 if (!madeTable)
1179                 {
1180                     bool isPlural = (ftype2 != -1);
1181                     gmx_fatal(FARGS, "Tabulated interaction of type '%s%s%s' with index %d cannot be used because no table file whose name matched '%s' was passed via the gmx mdrun -tableb command-line option.",
1182                               interaction_function[ftype1].longname,
1183                               isPlural ? "' or '" : "",
1184                               isPlural ? interaction_function[ftype2].longname : "",
1185                               i,
1186                               patternToFind.c_str());
1187                 }
1188             }
1189         }
1190         sfree(count);
1191     }
1192
1193     return tab;
1194 }
1195
1196 void forcerec_set_ranges(t_forcerec *fr,
1197                          int ncg_home, int ncg_force,
1198                          int natoms_force,
1199                          int natoms_force_constr, int natoms_f_novirsum)
1200 {
1201     fr->cg0 = 0;
1202     fr->hcg = ncg_home;
1203
1204     /* fr->ncg_force is unused in the standard code,
1205      * but it can be useful for modified code dealing with charge groups.
1206      */
1207     fr->ncg_force           = ncg_force;
1208     fr->natoms_force        = natoms_force;
1209     fr->natoms_force_constr = natoms_force_constr;
1210
1211     if (fr->natoms_force_constr > fr->nalloc_force)
1212     {
1213         fr->nalloc_force = over_alloc_dd(fr->natoms_force_constr);
1214     }
1215
1216     if (fr->haveDirectVirialContributions)
1217     {
1218         fr->forceBufferForDirectVirialContributions->resize(natoms_f_novirsum);
1219     }
1220 }
1221
1222 static real cutoff_inf(real cutoff)
1223 {
1224     if (cutoff == 0)
1225     {
1226         cutoff = GMX_CUTOFF_INF;
1227     }
1228
1229     return cutoff;
1230 }
1231
1232 /*! \brief Print Coulomb Ewald citations and set ewald coefficients */
1233 static void initCoulombEwaldParameters(FILE *fp, const t_inputrec *ir,
1234                                        bool systemHasNetCharge,
1235                                        interaction_const_t *ic)
1236 {
1237     if (!EEL_PME_EWALD(ir->coulombtype))
1238     {
1239         return;
1240     }
1241
1242     if (fp)
1243     {
1244         fprintf(fp, "Will do PME sum in reciprocal space for electrostatic interactions.\n");
1245
1246         if (ir->coulombtype == eelP3M_AD)
1247         {
1248             please_cite(fp, "Hockney1988");
1249             please_cite(fp, "Ballenegger2012");
1250         }
1251         else
1252         {
1253             please_cite(fp, "Essmann95a");
1254         }
1255
1256         if (ir->ewald_geometry == eewg3DC)
1257         {
1258             if (fp)
1259             {
1260                 fprintf(fp, "Using the Ewald3DC correction for systems with a slab geometry%s.\n",
1261                         systemHasNetCharge ? " and net charge" : "");
1262             }
1263             please_cite(fp, "In-Chul99a");
1264             if (systemHasNetCharge)
1265             {
1266                 please_cite(fp, "Ballenegger2009");
1267             }
1268         }
1269     }
1270
1271     ic->ewaldcoeff_q = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
1272     if (fp)
1273     {
1274         fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for Ewald\n",
1275                 1/ic->ewaldcoeff_q);
1276     }
1277
1278     if (ic->coulomb_modifier == eintmodPOTSHIFT)
1279     {
1280         GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero");
1281         ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb) / ic->rcoulomb;
1282     }
1283     else
1284     {
1285         ic->sh_ewald = 0;
1286     }
1287 }
1288
1289 /*! \brief Print Van der Waals Ewald citations and set ewald coefficients */
1290 static void initVdwEwaldParameters(FILE *fp, const t_inputrec *ir,
1291                                    interaction_const_t *ic)
1292 {
1293     if (!EVDW_PME(ir->vdwtype))
1294     {
1295         return;
1296     }
1297
1298     if (fp)
1299     {
1300         fprintf(fp, "Will do PME sum in reciprocal space for LJ dispersion interactions.\n");
1301         please_cite(fp, "Essmann95a");
1302     }
1303     ic->ewaldcoeff_lj = calc_ewaldcoeff_lj(ir->rvdw, ir->ewald_rtol_lj);
1304     if (fp)
1305     {
1306         fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for LJ Ewald\n",
1307                 1/ic->ewaldcoeff_lj);
1308     }
1309
1310     if (ic->vdw_modifier == eintmodPOTSHIFT)
1311     {
1312         real crc2       = gmx::square(ic->ewaldcoeff_lj*ic->rvdw);
1313         ic->sh_lj_ewald = (std::exp(-crc2)*(1 + crc2 + 0.5*crc2*crc2) - 1)/gmx::power6(ic->rvdw);
1314     }
1315     else
1316     {
1317         ic->sh_lj_ewald = 0;
1318     }
1319 }
1320
1321 gmx_bool uses_simple_tables(int                       cutoff_scheme,
1322                             const nonbonded_verlet_t *nbv)
1323 {
1324     gmx_bool bUsesSimpleTables = TRUE;
1325
1326     switch (cutoff_scheme)
1327     {
1328         case ecutsGROUP:
1329             bUsesSimpleTables = TRUE;
1330             break;
1331         case ecutsVERLET:
1332             GMX_RELEASE_ASSERT(nullptr != nbv, "A non-bonded verlet object is required with the Verlet cutoff-scheme");
1333             bUsesSimpleTables = nbv->pairlistIsSimple();
1334             break;
1335         default:
1336             gmx_incons("unimplemented");
1337     }
1338     return bUsesSimpleTables;
1339 }
1340
1341 static void init_ewald_f_table(interaction_const_t *ic,
1342                                real                 rtab)
1343 {
1344     real maxr;
1345
1346     /* Get the Ewald table spacing based on Coulomb and/or LJ
1347      * Ewald coefficients and rtol.
1348      */
1349     ic->tabq_scale = ewald_spline3_table_scale(ic);
1350
1351     if (ic->cutoff_scheme == ecutsVERLET)
1352     {
1353         maxr = ic->rcoulomb;
1354     }
1355     else
1356     {
1357         maxr = std::max(ic->rcoulomb, rtab);
1358     }
1359     ic->tabq_size  = static_cast<int>(maxr*ic->tabq_scale) + 2;
1360
1361     sfree_aligned(ic->tabq_coul_FDV0);
1362     sfree_aligned(ic->tabq_coul_F);
1363     sfree_aligned(ic->tabq_coul_V);
1364
1365     sfree_aligned(ic->tabq_vdw_FDV0);
1366     sfree_aligned(ic->tabq_vdw_F);
1367     sfree_aligned(ic->tabq_vdw_V);
1368
1369     if (EEL_PME_EWALD(ic->eeltype))
1370     {
1371         /* Create the original table data in FDV0 */
1372         snew_aligned(ic->tabq_coul_FDV0, ic->tabq_size*4, 32);
1373         snew_aligned(ic->tabq_coul_F, ic->tabq_size, 32);
1374         snew_aligned(ic->tabq_coul_V, ic->tabq_size, 32);
1375         table_spline3_fill_ewald_lr(ic->tabq_coul_F, ic->tabq_coul_V, ic->tabq_coul_FDV0,
1376                                     ic->tabq_size, 1/ic->tabq_scale, ic->ewaldcoeff_q, v_q_ewald_lr);
1377     }
1378
1379     if (EVDW_PME(ic->vdwtype))
1380     {
1381         snew_aligned(ic->tabq_vdw_FDV0, ic->tabq_size*4, 32);
1382         snew_aligned(ic->tabq_vdw_F, ic->tabq_size, 32);
1383         snew_aligned(ic->tabq_vdw_V, ic->tabq_size, 32);
1384         table_spline3_fill_ewald_lr(ic->tabq_vdw_F, ic->tabq_vdw_V, ic->tabq_vdw_FDV0,
1385                                     ic->tabq_size, 1/ic->tabq_scale, ic->ewaldcoeff_lj, v_lj_ewald_lr);
1386     }
1387 }
1388
1389 void init_interaction_const_tables(FILE                *fp,
1390                                    interaction_const_t *ic,
1391                                    real                 rtab)
1392 {
1393     if (EEL_PME_EWALD(ic->eeltype) || EVDW_PME(ic->vdwtype))
1394     {
1395         init_ewald_f_table(ic, rtab);
1396
1397         if (fp != nullptr)
1398         {
1399             fprintf(fp, "Initialized non-bonded Ewald correction tables, spacing: %.2e size: %d\n\n",
1400                     1/ic->tabq_scale, ic->tabq_size);
1401         }
1402     }
1403 }
1404
1405 static void clear_force_switch_constants(shift_consts_t *sc)
1406 {
1407     sc->c2   = 0;
1408     sc->c3   = 0;
1409     sc->cpot = 0;
1410 }
1411
1412 static void force_switch_constants(real p,
1413                                    real rsw, real rc,
1414                                    shift_consts_t *sc)
1415 {
1416     /* Here we determine the coefficient for shifting the force to zero
1417      * between distance rsw and the cut-off rc.
1418      * For a potential of r^-p, we have force p*r^-(p+1).
1419      * But to save flops we absorb p in the coefficient.
1420      * Thus we get:
1421      * force/p   = r^-(p+1) + c2*r^2 + c3*r^3
1422      * potential = r^-p + c2/3*r^3 + c3/4*r^4 + cpot
1423      */
1424     sc->c2   =  ((p + 1)*rsw - (p + 4)*rc)/(pow(rc, p + 2)*gmx::square(rc - rsw));
1425     sc->c3   = -((p + 1)*rsw - (p + 3)*rc)/(pow(rc, p + 2)*gmx::power3(rc - rsw));
1426     sc->cpot = -pow(rc, -p) + p*sc->c2/3*gmx::power3(rc - rsw) + p*sc->c3/4*gmx::power4(rc - rsw);
1427 }
1428
1429 static void potential_switch_constants(real rsw, real rc,
1430                                        switch_consts_t *sc)
1431 {
1432     /* The switch function is 1 at rsw and 0 at rc.
1433      * The derivative and second derivate are zero at both ends.
1434      * rsw        = max(r - r_switch, 0)
1435      * sw         = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5
1436      * dsw        = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4
1437      * force      = force*dsw - potential*sw
1438      * potential *= sw
1439      */
1440     sc->c3 = -10/gmx::power3(rc - rsw);
1441     sc->c4 =  15/gmx::power4(rc - rsw);
1442     sc->c5 =  -6/gmx::power5(rc - rsw);
1443 }
1444
1445 /*! \brief Construct interaction constants
1446  *
1447  * This data is used (particularly) by search and force code for
1448  * short-range interactions. Many of these are constant for the whole
1449  * simulation; some are constant only after PME tuning completes.
1450  */
1451 static void
1452 init_interaction_const(FILE                       *fp,
1453                        interaction_const_t       **interaction_const,
1454                        const t_inputrec           *ir,
1455                        const gmx_mtop_t           *mtop,
1456                        bool                        systemHasNetCharge)
1457 {
1458     interaction_const_t *ic;
1459
1460     snew(ic, 1);
1461
1462     ic->cutoff_scheme   = ir->cutoff_scheme;
1463
1464     /* Just allocate something so we can free it */
1465     snew_aligned(ic->tabq_coul_FDV0, 16, 32);
1466     snew_aligned(ic->tabq_coul_F, 16, 32);
1467     snew_aligned(ic->tabq_coul_V, 16, 32);
1468
1469     /* Lennard-Jones */
1470     ic->vdwtype         = ir->vdwtype;
1471     ic->vdw_modifier    = ir->vdw_modifier;
1472     ic->reppow          = mtop->ffparams.reppow;
1473     ic->rvdw            = cutoff_inf(ir->rvdw);
1474     ic->rvdw_switch     = ir->rvdw_switch;
1475     ic->ljpme_comb_rule = ir->ljpme_combination_rule;
1476     ic->useBuckingham   = (mtop->ffparams.functype[0] == F_BHAM);
1477     if (ic->useBuckingham)
1478     {
1479         ic->buckinghamBMax = calcBuckinghamBMax(fp, mtop);
1480     }
1481
1482     initVdwEwaldParameters(fp, ir, ic);
1483
1484     clear_force_switch_constants(&ic->dispersion_shift);
1485     clear_force_switch_constants(&ic->repulsion_shift);
1486
1487     switch (ic->vdw_modifier)
1488     {
1489         case eintmodPOTSHIFT:
1490             /* Only shift the potential, don't touch the force */
1491             ic->dispersion_shift.cpot = -1.0/gmx::power6(ic->rvdw);
1492             ic->repulsion_shift.cpot  = -1.0/gmx::power12(ic->rvdw);
1493             break;
1494         case eintmodFORCESWITCH:
1495             /* Switch the force, switch and shift the potential */
1496             force_switch_constants(6.0, ic->rvdw_switch, ic->rvdw,
1497                                    &ic->dispersion_shift);
1498             force_switch_constants(12.0, ic->rvdw_switch, ic->rvdw,
1499                                    &ic->repulsion_shift);
1500             break;
1501         case eintmodPOTSWITCH:
1502             /* Switch the potential and force */
1503             potential_switch_constants(ic->rvdw_switch, ic->rvdw,
1504                                        &ic->vdw_switch);
1505             break;
1506         case eintmodNONE:
1507         case eintmodEXACTCUTOFF:
1508             /* Nothing to do here */
1509             break;
1510         default:
1511             gmx_incons("unimplemented potential modifier");
1512     }
1513
1514     ic->sh_invrc6 = -ic->dispersion_shift.cpot;
1515
1516     /* Electrostatics */
1517     ic->eeltype          = ir->coulombtype;
1518     ic->coulomb_modifier = ir->coulomb_modifier;
1519     ic->rcoulomb         = cutoff_inf(ir->rcoulomb);
1520     ic->rcoulomb_switch  = ir->rcoulomb_switch;
1521     ic->epsilon_r        = ir->epsilon_r;
1522
1523     /* Set the Coulomb energy conversion factor */
1524     if (ic->epsilon_r != 0)
1525     {
1526         ic->epsfac = ONE_4PI_EPS0/ic->epsilon_r;
1527     }
1528     else
1529     {
1530         /* eps = 0 is infinite dieletric: no Coulomb interactions */
1531         ic->epsfac = 0;
1532     }
1533
1534     /* Reaction-field */
1535     if (EEL_RF(ic->eeltype))
1536     {
1537         ic->epsilon_rf = ir->epsilon_rf;
1538         /* Generalized reaction field parameters are updated every step */
1539         if (ic->eeltype != eelGRF)
1540         {
1541             calc_rffac(fp, ic->eeltype, ic->epsilon_r, ic->epsilon_rf,
1542                        ic->rcoulomb, 0, 0, nullptr,
1543                        &ic->k_rf, &ic->c_rf);
1544         }
1545
1546         if (ir->cutoff_scheme == ecutsGROUP && ic->eeltype == eelRF_ZERO)
1547         {
1548             /* grompp should have done this, but this scheme is obsolete */
1549             ic->coulomb_modifier = eintmodEXACTCUTOFF;
1550         }
1551     }
1552     else
1553     {
1554         /* For plain cut-off we might use the reaction-field kernels */
1555         ic->epsilon_rf = ic->epsilon_r;
1556         ic->k_rf       = 0;
1557         if (ir->coulomb_modifier == eintmodPOTSHIFT)
1558         {
1559             ic->c_rf   = 1/ic->rcoulomb;
1560         }
1561         else
1562         {
1563             ic->c_rf   = 0;
1564         }
1565     }
1566
1567     initCoulombEwaldParameters(fp, ir, systemHasNetCharge, ic);
1568
1569     if (fp != nullptr)
1570     {
1571         real dispersion_shift;
1572
1573         dispersion_shift = ic->dispersion_shift.cpot;
1574         if (EVDW_PME(ic->vdwtype))
1575         {
1576             dispersion_shift -= ic->sh_lj_ewald;
1577         }
1578         fprintf(fp, "Potential shift: LJ r^-12: %.3e r^-6: %.3e",
1579                 ic->repulsion_shift.cpot, dispersion_shift);
1580
1581         if (ic->eeltype == eelCUT)
1582         {
1583             fprintf(fp, ", Coulomb %.e", -ic->c_rf);
1584         }
1585         else if (EEL_PME(ic->eeltype))
1586         {
1587             fprintf(fp, ", Ewald %.3e", -ic->sh_ewald);
1588         }
1589         fprintf(fp, "\n");
1590     }
1591
1592     *interaction_const = ic;
1593 }
1594
1595 static void
1596 done_interaction_const(interaction_const_t *interaction_const)
1597 {
1598     sfree_aligned(interaction_const->tabq_coul_FDV0);
1599     sfree_aligned(interaction_const->tabq_coul_F);
1600     sfree_aligned(interaction_const->tabq_coul_V);
1601     sfree(interaction_const);
1602 }
1603
1604 void init_forcerec(FILE                             *fp,
1605                    const gmx::MDLogger              &mdlog,
1606                    t_forcerec                       *fr,
1607                    t_fcdata                         *fcd,
1608                    const t_inputrec                 *ir,
1609                    const gmx_mtop_t                 *mtop,
1610                    const t_commrec                  *cr,
1611                    matrix                            box,
1612                    const char                       *tabfn,
1613                    const char                       *tabpfn,
1614                    gmx::ArrayRef<const std::string>  tabbfnm,
1615                    const gmx_hw_info_t              &hardwareInfo,
1616                    const gmx_device_info_t          *deviceInfo,
1617                    const bool                        useGpuForBonded,
1618                    gmx_bool                          bNoSolvOpt,
1619                    real                              print_force)
1620 {
1621     int            m, negp_pp, negptable, egi, egj;
1622     real           rtab;
1623     char          *env;
1624     double         dbl;
1625     const t_block *cgs;
1626     gmx_bool       bGenericKernelOnly;
1627     gmx_bool       needGroupSchemeTables, bSomeNormalNbListsAreInUse;
1628     gmx_bool       bFEP_NonBonded;
1629     int           *nm_ind, egp_flags;
1630
1631
1632     /* By default we turn SIMD kernels on, but it might be turned off further down... */
1633     fr->use_simd_kernels = TRUE;
1634
1635     if (check_box(ir->ePBC, box))
1636     {
1637         gmx_fatal(FARGS, "%s", check_box(ir->ePBC, box));
1638     }
1639
1640     /* Test particle insertion ? */
1641     if (EI_TPI(ir->eI))
1642     {
1643         /* Set to the size of the molecule to be inserted (the last one) */
1644         /* Because of old style topologies, we have to use the last cg
1645          * instead of the last molecule type.
1646          */
1647         cgs       = &mtop->moltype[mtop->molblock.back().type].cgs;
1648         fr->n_tpi = cgs->index[cgs->nr] - cgs->index[cgs->nr-1];
1649         gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
1650         if (fr->n_tpi != molecules.block(molecules.numBlocks() - 1).size())
1651         {
1652             gmx_fatal(FARGS, "The molecule to insert can not consist of multiple charge groups.\nMake it a single charge group.");
1653         }
1654     }
1655     else
1656     {
1657         fr->n_tpi = 0;
1658     }
1659
1660     if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
1661     {
1662         gmx_fatal(FARGS, "%s electrostatics is no longer supported",
1663                   eel_names[ir->coulombtype]);
1664     }
1665
1666     if (ir->bAdress)
1667     {
1668         gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1669     }
1670     if (ir->useTwinRange)
1671     {
1672         gmx_fatal(FARGS, "Twin-range simulations are no longer supported");
1673     }
1674     /* Copy the user determined parameters */
1675     fr->userint1  = ir->userint1;
1676     fr->userint2  = ir->userint2;
1677     fr->userint3  = ir->userint3;
1678     fr->userint4  = ir->userint4;
1679     fr->userreal1 = ir->userreal1;
1680     fr->userreal2 = ir->userreal2;
1681     fr->userreal3 = ir->userreal3;
1682     fr->userreal4 = ir->userreal4;
1683
1684     /* Shell stuff */
1685     fr->fc_stepsize = ir->fc_stepsize;
1686
1687     /* Free energy */
1688     fr->efep        = ir->efep;
1689     fr->sc_alphavdw = ir->fepvals->sc_alpha;
1690     if (ir->fepvals->bScCoul)
1691     {
1692         fr->sc_alphacoul  = ir->fepvals->sc_alpha;
1693         fr->sc_sigma6_min = gmx::power6(ir->fepvals->sc_sigma_min);
1694     }
1695     else
1696     {
1697         fr->sc_alphacoul  = 0;
1698         fr->sc_sigma6_min = 0; /* only needed when bScCoul is on */
1699     }
1700     fr->sc_power      = ir->fepvals->sc_power;
1701     fr->sc_r_power    = ir->fepvals->sc_r_power;
1702     fr->sc_sigma6_def = gmx::power6(ir->fepvals->sc_sigma);
1703
1704     env = getenv("GMX_SCSIGMA_MIN");
1705     if (env != nullptr)
1706     {
1707         dbl = 0;
1708         sscanf(env, "%20lf", &dbl);
1709         fr->sc_sigma6_min = gmx::power6(dbl);
1710         if (fp)
1711         {
1712             fprintf(fp, "Setting the minimum soft core sigma to %g nm\n", dbl);
1713         }
1714     }
1715
1716     fr->bNonbonded = TRUE;
1717     if (getenv("GMX_NO_NONBONDED") != nullptr)
1718     {
1719         /* turn off non-bonded calculations */
1720         fr->bNonbonded = FALSE;
1721         GMX_LOG(mdlog.warning).asParagraph().appendText(
1722                 "Found environment variable GMX_NO_NONBONDED.\n"
1723                 "Disabling nonbonded calculations.");
1724     }
1725
1726     bGenericKernelOnly = FALSE;
1727
1728     /* We now check in the NS code whether a particular combination of interactions
1729      * can be used with water optimization, and disable it if that is not the case.
1730      */
1731
1732     if (getenv("GMX_NB_GENERIC") != nullptr)
1733     {
1734         if (fp != nullptr)
1735         {
1736             fprintf(fp,
1737                     "Found environment variable GMX_NB_GENERIC.\n"
1738                     "Disabling all interaction-specific nonbonded kernels, will only\n"
1739                     "use the slow generic ones in src/gmxlib/nonbonded/nb_generic.c\n\n");
1740         }
1741         bGenericKernelOnly = TRUE;
1742     }
1743
1744     if (bGenericKernelOnly)
1745     {
1746         bNoSolvOpt         = TRUE;
1747     }
1748
1749     if ( (getenv("GMX_DISABLE_SIMD_KERNELS") != nullptr) || (getenv("GMX_NOOPTIMIZEDKERNELS") != nullptr) )
1750     {
1751         fr->use_simd_kernels = FALSE;
1752         if (fp != nullptr)
1753         {
1754             fprintf(fp,
1755                     "\nFound environment variable GMX_DISABLE_SIMD_KERNELS.\n"
1756                     "Disabling the usage of any SIMD-specific non-bonded & bonded kernel routines\n"
1757                     "(e.g. SSE2/SSE4.1/AVX).\n\n");
1758         }
1759     }
1760
1761     fr->bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
1762
1763     /* Neighbour searching stuff */
1764     fr->cutoff_scheme = ir->cutoff_scheme;
1765     fr->bGrid         = (ir->ns_type == ensGRID);
1766     fr->ePBC          = ir->ePBC;
1767
1768     if (fr->cutoff_scheme == ecutsGROUP)
1769     {
1770         const char *note = "NOTE: This file uses the deprecated 'group' cutoff_scheme. This will be\n"
1771             "removed in a future release when 'verlet' supports all interaction forms.\n";
1772
1773         if (MASTER(cr))
1774         {
1775             fprintf(stderr, "\n%s\n", note);
1776         }
1777         if (fp != nullptr)
1778         {
1779             fprintf(fp, "\n%s\n", note);
1780         }
1781     }
1782
1783     /* Determine if we will do PBC for distances in bonded interactions */
1784     if (fr->ePBC == epbcNONE)
1785     {
1786         fr->bMolPBC = FALSE;
1787     }
1788     else
1789     {
1790         if (!DOMAINDECOMP(cr))
1791         {
1792             gmx_bool bSHAKE;
1793
1794             bSHAKE = (ir->eConstrAlg == econtSHAKE &&
1795                       (gmx_mtop_ftype_count(mtop, F_CONSTR) > 0 ||
1796                        gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0));
1797
1798             /* The group cut-off scheme and SHAKE assume charge groups
1799              * are whole, but not using molpbc is faster in most cases.
1800              * With intermolecular interactions we need PBC for calculating
1801              * distances between atoms in different molecules.
1802              */
1803             if ((fr->cutoff_scheme == ecutsGROUP || bSHAKE) &&
1804                 !mtop->bIntermolecularInteractions)
1805             {
1806                 fr->bMolPBC = ir->bPeriodicMols;
1807
1808                 if (bSHAKE && fr->bMolPBC)
1809                 {
1810                     gmx_fatal(FARGS, "SHAKE is not supported with periodic molecules");
1811                 }
1812             }
1813             else
1814             {
1815                 /* Not making molecules whole is faster in most cases,
1816                  * but With orientation restraints we need whole molecules.
1817                  */
1818                 fr->bMolPBC = (fcd->orires.nr == 0);
1819
1820                 if (getenv("GMX_USE_GRAPH") != nullptr)
1821                 {
1822                     fr->bMolPBC = FALSE;
1823                     if (fp)
1824                     {
1825                         GMX_LOG(mdlog.warning).asParagraph().appendText("GMX_USE_GRAPH is set, using the graph for bonded interactions");
1826                     }
1827
1828                     if (mtop->bIntermolecularInteractions)
1829                     {
1830                         GMX_LOG(mdlog.warning).asParagraph().appendText("WARNING: Molecules linked by intermolecular interactions have to reside in the same periodic image, otherwise artifacts will occur!");
1831                     }
1832                 }
1833
1834                 GMX_RELEASE_ASSERT(fr->bMolPBC || !mtop->bIntermolecularInteractions, "We need to use PBC within molecules with inter-molecular interactions");
1835
1836                 if (bSHAKE && fr->bMolPBC)
1837                 {
1838                     gmx_fatal(FARGS, "SHAKE is not properly supported with intermolecular interactions. For short simulations where linked molecules remain in the same periodic image, the environment variable GMX_USE_GRAPH can be used to override this check.\n");
1839                 }
1840             }
1841         }
1842         else
1843         {
1844             fr->bMolPBC = dd_bonded_molpbc(cr->dd, fr->ePBC);
1845         }
1846     }
1847
1848     fr->rc_scaling = ir->refcoord_scaling;
1849     copy_rvec(ir->posres_com, fr->posres_com);
1850     copy_rvec(ir->posres_comB, fr->posres_comB);
1851     fr->rlist                    = cutoff_inf(ir->rlist);
1852     fr->ljpme_combination_rule   = ir->ljpme_combination_rule;
1853
1854     /* This now calculates sum for q and c6*/
1855     bool systemHasNetCharge = set_chargesum(fp, fr, mtop);
1856
1857     /* fr->ic is used both by verlet and group kernels (to some extent) now */
1858     init_interaction_const(fp, &fr->ic, ir, mtop, systemHasNetCharge);
1859     init_interaction_const_tables(fp, fr->ic, ir->rlist + ir->tabext);
1860
1861     const interaction_const_t *ic = fr->ic;
1862
1863     /* TODO: Replace this Ewald table or move it into interaction_const_t */
1864     if (ir->coulombtype == eelEWALD)
1865     {
1866         init_ewald_tab(&(fr->ewald_table), ir, fp);
1867     }
1868
1869     /* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */
1870     switch (ic->eeltype)
1871     {
1872         case eelCUT:
1873             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_COULOMB;
1874             break;
1875
1876         case eelRF:
1877         case eelGRF:
1878             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_REACTIONFIELD;
1879             break;
1880
1881         case eelRF_ZERO:
1882             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_REACTIONFIELD;
1883             GMX_RELEASE_ASSERT(ic->coulomb_modifier == eintmodEXACTCUTOFF, "With the group scheme RF-zero needs the exact cut-off modifier");
1884             break;
1885
1886         case eelSWITCH:
1887         case eelSHIFT:
1888         case eelUSER:
1889         case eelENCADSHIFT:
1890         case eelPMESWITCH:
1891         case eelPMEUSER:
1892         case eelPMEUSERSWITCH:
1893             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_CUBICSPLINETABLE;
1894             break;
1895
1896         case eelPME:
1897         case eelP3M_AD:
1898         case eelEWALD:
1899             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_EWALD;
1900             break;
1901
1902         default:
1903             gmx_fatal(FARGS, "Unsupported electrostatic interaction: %s", eel_names[ic->eeltype]);
1904     }
1905     fr->nbkernel_elec_modifier = ic->coulomb_modifier;
1906
1907     /* Vdw: Translate from mdp settings to kernel format */
1908     switch (ic->vdwtype)
1909     {
1910         case evdwCUT:
1911             if (fr->bBHAM)
1912             {
1913                 fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_BUCKINGHAM;
1914             }
1915             else
1916             {
1917                 fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_LENNARDJONES;
1918             }
1919             break;
1920         case evdwPME:
1921             fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_LJEWALD;
1922             break;
1923
1924         case evdwSWITCH:
1925         case evdwSHIFT:
1926         case evdwUSER:
1927         case evdwENCADSHIFT:
1928             fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_CUBICSPLINETABLE;
1929             break;
1930
1931         default:
1932             gmx_fatal(FARGS, "Unsupported vdw interaction: %s", evdw_names[ic->vdwtype]);
1933     }
1934     fr->nbkernel_vdw_modifier = ic->vdw_modifier;
1935
1936     if (ir->cutoff_scheme == ecutsGROUP)
1937     {
1938         fr->bvdwtab    = ((ic->vdwtype != evdwCUT || !gmx_within_tol(ic->reppow, 12.0, 10*GMX_DOUBLE_EPS))
1939                           && !EVDW_PME(ic->vdwtype));
1940         /* We have special kernels for standard Ewald and PME, but the pme-switch ones are tabulated above */
1941         fr->bcoultab   = !(ic->eeltype == eelCUT ||
1942                            ic->eeltype == eelEWALD ||
1943                            ic->eeltype == eelPME ||
1944                            ic->eeltype == eelP3M_AD ||
1945                            ic->eeltype == eelRF ||
1946                            ic->eeltype == eelRF_ZERO);
1947
1948         /* If the user absolutely wants different switch/shift settings for coul/vdw, it is likely
1949          * going to be faster to tabulate the interaction than calling the generic kernel.
1950          * However, if generic kernels have been requested we keep things analytically.
1951          */
1952         if (fr->nbkernel_elec_modifier == eintmodPOTSWITCH &&
1953             fr->nbkernel_vdw_modifier == eintmodPOTSWITCH &&
1954             !bGenericKernelOnly)
1955         {
1956             if ((ic->rcoulomb_switch != ic->rvdw_switch) || (ic->rcoulomb != ic->rvdw))
1957             {
1958                 fr->bcoultab = TRUE;
1959                 /* Once we tabulate electrostatics, we can use the switch function for LJ,
1960                  * which would otherwise need two tables.
1961                  */
1962             }
1963         }
1964         else if ((fr->nbkernel_elec_modifier == eintmodPOTSHIFT && fr->nbkernel_vdw_modifier == eintmodPOTSHIFT) ||
1965                  ((fr->nbkernel_elec_interaction == GMX_NBKERNEL_ELEC_REACTIONFIELD &&
1966                    fr->nbkernel_elec_modifier == eintmodEXACTCUTOFF &&
1967                    (fr->nbkernel_vdw_modifier == eintmodPOTSWITCH || fr->nbkernel_vdw_modifier == eintmodPOTSHIFT))))
1968         {
1969             if ((ic->rcoulomb != ic->rvdw) && (!bGenericKernelOnly))
1970             {
1971                 fr->bcoultab = TRUE;
1972             }
1973         }
1974
1975         if (fr->nbkernel_elec_modifier == eintmodFORCESWITCH)
1976         {
1977             fr->bcoultab = TRUE;
1978         }
1979         if (fr->nbkernel_vdw_modifier == eintmodFORCESWITCH)
1980         {
1981             fr->bvdwtab = TRUE;
1982         }
1983
1984         if (getenv("GMX_REQUIRE_TABLES"))
1985         {
1986             fr->bvdwtab  = TRUE;
1987             fr->bcoultab = TRUE;
1988         }
1989
1990         if (fp)
1991         {
1992             fprintf(fp, "Table routines are used for coulomb: %s\n",
1993                     gmx::boolToString(fr->bcoultab));
1994             fprintf(fp, "Table routines are used for vdw:     %s\n",
1995                     gmx::boolToString(fr->bvdwtab));
1996         }
1997
1998         if (fr->bvdwtab)
1999         {
2000             fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_CUBICSPLINETABLE;
2001             fr->nbkernel_vdw_modifier    = eintmodNONE;
2002         }
2003         if (fr->bcoultab)
2004         {
2005             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_CUBICSPLINETABLE;
2006             fr->nbkernel_elec_modifier    = eintmodNONE;
2007         }
2008     }
2009
2010     if (ir->cutoff_scheme == ecutsVERLET)
2011     {
2012         if (!gmx_within_tol(ic->reppow, 12.0, 10*GMX_DOUBLE_EPS))
2013         {
2014             gmx_fatal(FARGS, "Cut-off scheme %s only supports LJ repulsion power 12", ecutscheme_names[ir->cutoff_scheme]);
2015         }
2016         /* Older tpr files can contain Coulomb user tables with the Verlet cutoff-scheme,
2017          * while mdrun does not (and never did) support this.
2018          */
2019         if (EEL_USER(fr->ic->eeltype))
2020         {
2021             gmx_fatal(FARGS, "Combination of %s and cutoff scheme %s is not supported",
2022                       eel_names[ir->coulombtype], ecutscheme_names[ir->cutoff_scheme]);
2023         }
2024
2025         fr->bvdwtab  = FALSE;
2026         fr->bcoultab = FALSE;
2027     }
2028
2029     /* 1-4 interaction electrostatics */
2030     fr->fudgeQQ = mtop->ffparams.fudgeQQ;
2031
2032     /* Parameters for generalized RF */
2033     fr->zsquare = 0.0;
2034     fr->temp    = 0.0;
2035
2036     if (ic->eeltype == eelGRF)
2037     {
2038         init_generalized_rf(fp, mtop, ir, fr);
2039     }
2040
2041     fr->haveDirectVirialContributions =
2042         (EEL_FULL(ic->eeltype) || EVDW_PME(ic->vdwtype) ||
2043          fr->forceProviders->hasForceProvider() ||
2044          gmx_mtop_ftype_count(mtop, F_POSRES) > 0 ||
2045          gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0 ||
2046          ir->nwall > 0 ||
2047          ir->bPull ||
2048          ir->bRot ||
2049          ir->bIMD);
2050
2051     if (fr->haveDirectVirialContributions)
2052     {
2053         fr->forceBufferForDirectVirialContributions = new std::vector<gmx::RVec>;
2054     }
2055
2056     if (fr->cutoff_scheme == ecutsGROUP &&
2057         ncg_mtop(mtop) > fr->cg_nalloc && !DOMAINDECOMP(cr))
2058     {
2059         /* Count the total number of charge groups */
2060         fr->cg_nalloc = ncg_mtop(mtop);
2061         srenew(fr->cg_cm, fr->cg_nalloc);
2062     }
2063     if (fr->shift_vec == nullptr)
2064     {
2065         snew(fr->shift_vec, SHIFTS);
2066     }
2067
2068     if (fr->fshift == nullptr)
2069     {
2070         snew(fr->fshift, SHIFTS);
2071     }
2072
2073     if (fr->nbfp == nullptr)
2074     {
2075         fr->ntype = mtop->ffparams.atnr;
2076         fr->nbfp  = mk_nbfp(&mtop->ffparams, fr->bBHAM);
2077         if (EVDW_PME(ic->vdwtype))
2078         {
2079             fr->ljpme_c6grid  = make_ljpme_c6grid(&mtop->ffparams, fr);
2080         }
2081     }
2082
2083     /* Copy the energy group exclusions */
2084     fr->egp_flags = ir->opts.egp_flags;
2085
2086     /* Van der Waals stuff */
2087     if ((ic->vdwtype != evdwCUT) && (ic->vdwtype != evdwUSER) && !fr->bBHAM)
2088     {
2089         if (ic->rvdw_switch >= ic->rvdw)
2090         {
2091             gmx_fatal(FARGS, "rvdw_switch (%f) must be < rvdw (%f)",
2092                       ic->rvdw_switch, ic->rvdw);
2093         }
2094         if (fp)
2095         {
2096             fprintf(fp, "Using %s Lennard-Jones, switch between %g and %g nm\n",
2097                     (ic->eeltype == eelSWITCH) ? "switched" : "shifted",
2098                     ic->rvdw_switch, ic->rvdw);
2099         }
2100     }
2101
2102     if (fr->bBHAM && EVDW_PME(ic->vdwtype))
2103     {
2104         gmx_fatal(FARGS, "LJ PME not supported with Buckingham");
2105     }
2106
2107     if (fr->bBHAM && (ic->vdwtype == evdwSHIFT || ic->vdwtype == evdwSWITCH))
2108     {
2109         gmx_fatal(FARGS, "Switch/shift interaction not supported with Buckingham");
2110     }
2111
2112     if (fr->bBHAM && fr->cutoff_scheme == ecutsVERLET)
2113     {
2114         gmx_fatal(FARGS, "Verlet cutoff-scheme is not supported with Buckingham");
2115     }
2116
2117     if (fp && fr->cutoff_scheme == ecutsGROUP)
2118     {
2119         fprintf(fp, "Cut-off's:   NS: %g   Coulomb: %g   %s: %g\n",
2120                 fr->rlist, ic->rcoulomb, fr->bBHAM ? "BHAM" : "LJ", ic->rvdw);
2121     }
2122
2123     fr->eDispCorr = ir->eDispCorr;
2124     fr->numAtomsForDispersionCorrection = mtop->natoms;
2125     if (ir->eDispCorr != edispcNO)
2126     {
2127         set_avcsixtwelve(fp, fr, mtop);
2128     }
2129
2130     if (ir->implicit_solvent)
2131     {
2132         gmx_fatal(FARGS, "Implict solvation is no longer supported.");
2133     }
2134
2135     /* Construct tables for the group scheme. A little unnecessary to
2136      * make both vdw and coul tables sometimes, but what the
2137      * heck. Note that both cutoff schemes construct Ewald tables in
2138      * init_interaction_const_tables. */
2139     needGroupSchemeTables = (ir->cutoff_scheme == ecutsGROUP &&
2140                              (fr->bcoultab || fr->bvdwtab));
2141
2142     negp_pp   = ir->opts.ngener - ir->nwall;
2143     negptable = 0;
2144     if (!needGroupSchemeTables)
2145     {
2146         bSomeNormalNbListsAreInUse = TRUE;
2147         fr->nnblists               = 1;
2148     }
2149     else
2150     {
2151         bSomeNormalNbListsAreInUse = FALSE;
2152         for (egi = 0; egi < negp_pp; egi++)
2153         {
2154             for (egj = egi; egj < negp_pp; egj++)
2155             {
2156                 egp_flags = ir->opts.egp_flags[GID(egi, egj, ir->opts.ngener)];
2157                 if (!(egp_flags & EGP_EXCL))
2158                 {
2159                     if (egp_flags & EGP_TABLE)
2160                     {
2161                         negptable++;
2162                     }
2163                     else
2164                     {
2165                         bSomeNormalNbListsAreInUse = TRUE;
2166                     }
2167                 }
2168             }
2169         }
2170         if (bSomeNormalNbListsAreInUse)
2171         {
2172             fr->nnblists = negptable + 1;
2173         }
2174         else
2175         {
2176             fr->nnblists = negptable;
2177         }
2178         if (fr->nnblists > 1)
2179         {
2180             snew(fr->gid2nblists, ir->opts.ngener*ir->opts.ngener);
2181         }
2182     }
2183
2184     snew(fr->nblists, fr->nnblists);
2185
2186     /* This code automatically gives table length tabext without cut-off's,
2187      * in that case grompp should already have checked that we do not need
2188      * normal tables and we only generate tables for 1-4 interactions.
2189      */
2190     rtab = ir->rlist + ir->tabext;
2191
2192     if (needGroupSchemeTables)
2193     {
2194         /* make tables for ordinary interactions */
2195         if (bSomeNormalNbListsAreInUse)
2196         {
2197             make_nbf_tables(fp, ic, rtab, tabfn, nullptr, nullptr, &fr->nblists[0]);
2198             m = 1;
2199         }
2200         else
2201         {
2202             m = 0;
2203         }
2204         if (negptable > 0)
2205         {
2206             /* Read the special tables for certain energy group pairs */
2207             nm_ind = mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].nm_ind;
2208             for (egi = 0; egi < negp_pp; egi++)
2209             {
2210                 for (egj = egi; egj < negp_pp; egj++)
2211                 {
2212                     egp_flags = ir->opts.egp_flags[GID(egi, egj, ir->opts.ngener)];
2213                     if ((egp_flags & EGP_TABLE) && !(egp_flags & EGP_EXCL))
2214                     {
2215                         if (fr->nnblists > 1)
2216                         {
2217                             fr->gid2nblists[GID(egi, egj, ir->opts.ngener)] = m;
2218                         }
2219                         /* Read the table file with the two energy groups names appended */
2220                         make_nbf_tables(fp, ic, rtab, tabfn,
2221                                         *mtop->groups.groupNames[nm_ind[egi]],
2222                                         *mtop->groups.groupNames[nm_ind[egj]],
2223                                         &fr->nblists[m]);
2224                         m++;
2225                     }
2226                     else if (fr->nnblists > 1)
2227                     {
2228                         fr->gid2nblists[GID(egi, egj, ir->opts.ngener)] = 0;
2229                     }
2230                 }
2231             }
2232         }
2233     }
2234
2235     /* Tables might not be used for the potential modifier
2236      * interactions per se, but we still need them to evaluate
2237      * switch/shift dispersion corrections in this case. */
2238     if (fr->eDispCorr != edispcNO)
2239     {
2240         fr->dispersionCorrectionTable = makeDispersionCorrectionTable(fp, ic, rtab, tabfn);
2241     }
2242
2243     /* We want to use unmodified tables for 1-4 coulombic
2244      * interactions, so we must in general have an extra set of
2245      * tables. */
2246     if (gmx_mtop_ftype_count(mtop, F_LJ14) > 0 ||
2247         gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0 ||
2248         gmx_mtop_ftype_count(mtop, F_LJC_PAIRS_NB) > 0)
2249     {
2250         fr->pairsTable = make_tables(fp, ic, tabpfn, rtab,
2251                                      GMX_MAKETABLES_14ONLY);
2252     }
2253
2254     /* Wall stuff */
2255     fr->nwall = ir->nwall;
2256     if (ir->nwall && ir->wall_type == ewtTABLE)
2257     {
2258         make_wall_tables(fp, ir, tabfn, &mtop->groups, fr);
2259     }
2260
2261     if (fcd && !tabbfnm.empty())
2262     {
2263         // Need to catch std::bad_alloc
2264         // TODO Don't need to catch this here, when merging with master branch
2265         try
2266         {
2267             fcd->bondtab  = make_bonded_tables(fp,
2268                                                F_TABBONDS, F_TABBONDSNC,
2269                                                mtop, tabbfnm, "b");
2270             fcd->angletab = make_bonded_tables(fp,
2271                                                F_TABANGLES, -1,
2272                                                mtop, tabbfnm, "a");
2273             fcd->dihtab   = make_bonded_tables(fp,
2274                                                F_TABDIHS, -1,
2275                                                mtop, tabbfnm, "d");
2276         }
2277         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2278     }
2279     else
2280     {
2281         if (debug)
2282         {
2283             fprintf(debug, "No fcdata or table file name passed, can not read table, can not do bonded interactions\n");
2284         }
2285     }
2286
2287     // QM/MM initialization if requested
2288     fr->bQMMM = ir->bQMMM;
2289     if (fr->bQMMM)
2290     {
2291         // Initialize QM/MM if supported
2292         if (GMX_QMMM)
2293         {
2294             GMX_LOG(mdlog.info).asParagraph().
2295                 appendText("Large parts of the QM/MM support is deprecated, and may be removed in a future "
2296                            "version. Please get in touch with the developers if you find the support useful, "
2297                            "as help is needed if the functionality is to continue to be available.");
2298             fr->qr = mk_QMMMrec();
2299             init_QMMMrec(cr, mtop, ir, fr);
2300         }
2301         else
2302         {
2303             gmx_incons("QM/MM was requested, but is only available when GROMACS "
2304                        "is configured with QM/MM support");
2305         }
2306     }
2307
2308     /* Set all the static charge group info */
2309     fr->cginfo_mb = init_cginfo_mb(fp, mtop, fr, bNoSolvOpt,
2310                                    &bFEP_NonBonded,
2311                                    &fr->bExcl_IntraCGAll_InterCGNone);
2312     if (DOMAINDECOMP(cr))
2313     {
2314         fr->cginfo = nullptr;
2315     }
2316     else
2317     {
2318         fr->cginfo = cginfo_expand(mtop->molblock.size(), fr->cginfo_mb);
2319     }
2320
2321     if (!DOMAINDECOMP(cr))
2322     {
2323         forcerec_set_ranges(fr, ncg_mtop(mtop), ncg_mtop(mtop),
2324                             mtop->natoms, mtop->natoms, mtop->natoms);
2325     }
2326
2327     fr->print_force = print_force;
2328
2329
2330     /* Initialize neighbor search */
2331     snew(fr->ns, 1);
2332     init_ns(fp, cr, fr->ns, fr, mtop);
2333
2334     /* Initialize the thread working data for bonded interactions */
2335     init_bonded_threading(fp, mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].nr,
2336                           &fr->bondedThreading);
2337
2338     fr->nthread_ewc = gmx_omp_nthreads_get(emntBonded);
2339     snew(fr->ewc_t, fr->nthread_ewc);
2340
2341     if (fr->cutoff_scheme == ecutsVERLET)
2342     {
2343         // We checked the cut-offs in grompp, but double-check here.
2344         // We have PME+LJcutoff kernels for rcoulomb>rvdw.
2345         if (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == eelCUT)
2346         {
2347             GMX_RELEASE_ASSERT(ir->rcoulomb >= ir->rvdw, "With Verlet lists and PME we should have rcoulomb>=rvdw");
2348         }
2349         else
2350         {
2351             GMX_RELEASE_ASSERT(ir->rcoulomb == ir->rvdw, "With Verlet lists and no PME rcoulomb and rvdw should be identical");
2352         }
2353
2354         fr->nbv = Nbnxm::init_nb_verlet(mdlog, bFEP_NonBonded, ir, fr,
2355                                         cr, hardwareInfo, deviceInfo,
2356                                         mtop, box);
2357
2358         if (useGpuForBonded)
2359         {
2360             auto stream = DOMAINDECOMP(cr) ?
2361                 Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, Nbnxm::InteractionLocality::NonLocal) :
2362                 Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, Nbnxm::InteractionLocality::Local);
2363             // TODO the heap allocation is only needed while
2364             // t_forcerec lacks a constructor.
2365             fr->gpuBonded = new gmx::GpuBonded(mtop->ffparams,
2366                                                stream);
2367         }
2368     }
2369
2370     if (fp != nullptr)
2371     {
2372         /* Here we switch from using mdlog, which prints the newline before
2373          * the paragraph, to our old fprintf logging, which prints the newline
2374          * after the paragraph, so we should add a newline here.
2375          */
2376         fprintf(fp, "\n");
2377     }
2378
2379     if (ir->eDispCorr != edispcNO)
2380     {
2381         calc_enervirdiff(fp, ir->eDispCorr, fr);
2382     }
2383 }
2384
2385 /* Frees GPU memory and sets a tMPI node barrier.
2386  *
2387  * Note that this function needs to be called even if GPUs are not used
2388  * in this run because the PME ranks have no knowledge of whether GPUs
2389  * are used or not, but all ranks need to enter the barrier below.
2390  * \todo Remove physical node barrier from this function after making sure
2391  * that it's not needed anymore (with a shared GPU run).
2392  */
2393 void free_gpu_resources(t_forcerec                          *fr,
2394                         const gmx::PhysicalNodeCommunicator &physicalNodeCommunicator)
2395 {
2396     bool isPPrankUsingGPU = (fr != nullptr) && (fr->nbv != nullptr) && fr->nbv->useGpu();
2397
2398     /* stop the GPU profiler (only CUDA) */
2399     stopGpuProfiler();
2400
2401     if (isPPrankUsingGPU)
2402     {
2403         /* Free data in GPU memory and pinned memory before destroying the GPU context */
2404         fr->nbv.reset();
2405
2406         delete fr->gpuBonded;
2407         fr->gpuBonded = nullptr;
2408     }
2409
2410     /* With tMPI we need to wait for all ranks to finish deallocation before
2411      * destroying the CUDA context in free_gpu() as some tMPI ranks may be sharing
2412      * GPU and context.
2413      *
2414      * This is not a concern in OpenCL where we use one context per rank which
2415      * is freed in nbnxn_gpu_free().
2416      *
2417      * Note: it is safe to not call the barrier on the ranks which do not use GPU,
2418      * but it is easier and more futureproof to call it on the whole node.
2419      */
2420     if (GMX_THREAD_MPI)
2421     {
2422         physicalNodeCommunicator.barrier();
2423     }
2424 }
2425
2426 void done_forcerec(t_forcerec *fr, int numMolBlocks, int numEnergyGroups)
2427 {
2428     if (fr == nullptr)
2429     {
2430         // PME-only ranks don't have a forcerec
2431         return;
2432     }
2433     // cginfo is dynamically allocated if no domain decomposition
2434     if (fr->cginfo != nullptr)
2435     {
2436         sfree(fr->cginfo);
2437     }
2438     done_cginfo_mb(fr->cginfo_mb, numMolBlocks);
2439     sfree(fr->nbfp);
2440     done_interaction_const(fr->ic);
2441     sfree(fr->shift_vec);
2442     sfree(fr->fshift);
2443     sfree(fr->nblists);
2444     done_ns(fr->ns, numEnergyGroups);
2445     sfree(fr->ewc_t);
2446     tear_down_bonded_threading(fr->bondedThreading);
2447     GMX_RELEASE_ASSERT(fr->gpuBonded == nullptr, "Should have been deleted earlier, when used");
2448     fr->bondedThreading = nullptr;
2449     delete fr;
2450 }