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