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-2019,2020,2021, 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/ewald/pme_pp_comm_gpu.h"
  57 #include "gromacs/fileio/filetypes.h"
  58 #include "gromacs/gmxlib/network.h"
  59 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  60 #include "gromacs/gpu_utils/gpu_utils.h"
  61 #include "gromacs/hardware/hw_info.h"
  62 #include "gromacs/listed_forces/gpubonded.h"
  63 #include "gromacs/listed_forces/listed_forces.h"
  64 #include "gromacs/listed_forces/pairs.h"
  65 #include "gromacs/math/functions.h"
  66 #include "gromacs/math/units.h"
  67 #include "gromacs/math/utilities.h"
  68 #include "gromacs/math/vec.h"
  69 #include "gromacs/mdlib/dispersioncorrection.h"
  70 #include "gromacs/mdlib/force.h"
  71 #include "gromacs/mdlib/forcerec_threading.h"
  72 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  73 #include "gromacs/mdlib/md_support.h"
  74 #include "gromacs/mdlib/rf_util.h"
  75 #include "gromacs/mdlib/wall.h"
  76 #include "gromacs/mdlib/wholemoleculetransform.h"
  77 #include "gromacs/mdtypes/commrec.h"
  78 #include "gromacs/mdtypes/fcdata.h"
  79 #include "gromacs/mdtypes/forcerec.h"
  80 #include "gromacs/mdtypes/group.h"
  81 #include "gromacs/mdtypes/iforceprovider.h"
  82 #include "gromacs/mdtypes/inputrec.h"
  83 #include "gromacs/mdtypes/interaction_const.h"
  84 #include "gromacs/mdtypes/md_enums.h"
  85 #include "gromacs/mdtypes/multipletimestepping.h"
  86 #include "gromacs/mdtypes/nblist.h"
  87 #include "gromacs/nbnxm/nbnxm.h"
  88 #include "gromacs/pbcutil/ishift.h"
  89 #include "gromacs/pbcutil/pbc.h"
  90 #include "gromacs/tables/forcetable.h"
  91 #include "gromacs/topology/mtop_util.h"
  92 #include "gromacs/trajectory/trajectoryframe.h"
  93 #include "gromacs/utility/cstringutil.h"
  94 #include "gromacs/utility/exceptions.h"
  95 #include "gromacs/utility/fatalerror.h"
  96 #include "gromacs/utility/gmxassert.h"
  97 #include "gromacs/utility/logger.h"
  98 #include "gromacs/utility/physicalnodecommunicator.h"
  99 #include "gromacs/utility/pleasecite.h"
 100 #include "gromacs/utility/smalloc.h"
 101 #include "gromacs/utility/strconvert.h"
 102
 103 #include "gpuforcereduction.h"
 104
 105 ForceHelperBuffers::ForceHelperBuffers(bool haveDirectVirialContributions) :
 106     haveDirectVirialContributions_(haveDirectVirialContributions)
 107 {
 108     shiftForces_.resize(SHIFTS);
 109 }
 110
 111 void ForceHelperBuffers::resize(int numAtoms)
 112 {
 113     if (haveDirectVirialContributions_)
 114     {
 115         forceBufferForDirectVirialContributions_.resize(numAtoms);
 116     }
 117 }
 118
 119 std::vector<real> makeNonBondedParameterLists(const gmx_ffparams_t& forceFieldParams, bool useBuckinghamPotential)
 120 {
 121     std::vector<real> nbfp;
 122     int               atnr;
 123
 124     atnr = forceFieldParams.atnr;
 125     if (useBuckinghamPotential)
 126     {
 127         nbfp.resize(3 * atnr * atnr);
 128         int k = 0;
 129         for (int i = 0; (i < atnr); i++)
 130         {
 131             for (int j = 0; (j < atnr); j++, k++)
 132             {
 133                 BHAMA(nbfp, atnr, i, j) = forceFieldParams.iparams[k].bham.a;
 134                 BHAMB(nbfp, atnr, i, j) = forceFieldParams.iparams[k].bham.b;
 135                 /* nbfp now includes the 6.0 derivative prefactor */
 136                 BHAMC(nbfp, atnr, i, j) = forceFieldParams.iparams[k].bham.c * 6.0;
 137             }
 138         }
 139     }
 140     else
 141     {
 142         nbfp.resize(2 * atnr * atnr);
 143         int k = 0;
 144         for (int i = 0; (i < atnr); i++)
 145         {
 146             for (int j = 0; (j < atnr); j++, k++)
 147             {
 148                 /* nbfp now includes the 6.0/12.0 derivative prefactors */
 149                 C6(nbfp, atnr, i, j)  = forceFieldParams.iparams[k].lj.c6 * 6.0;
 150                 C12(nbfp, atnr, i, j) = forceFieldParams.iparams[k].lj.c12 * 12.0;
 151             }
 152         }
 153     }
 154
 155     return nbfp;
 156 }
 157
 158 std::vector<real> makeLJPmeC6GridCorrectionParameters(const gmx_ffparams_t& forceFieldParams,
 159                                                       const t_forcerec&     forceRec)
 160 {
 161     int  i, j, k, atnr;
 162     real c6, c6i, c6j, c12i, c12j, epsi, epsj, sigmai, sigmaj;
 163
 164     /* For LJ-PME simulations, we correct the energies with the reciprocal space
 165      * inside of the cut-off. To do this the non-bonded kernels needs to have
 166      * access to the C6-values used on the reciprocal grid in pme.c
 167      */
 168
 169     atnr = forceFieldParams.atnr;
 170     std::vector<real> grid(2 * atnr * atnr, 0.0);
 171     for (i = k = 0; (i < atnr); i++)
 172     {
 173         for (j = 0; (j < atnr); j++, k++)
 174         {
 175             c6i  = forceFieldParams.iparams[i * (atnr + 1)].lj.c6;
 176             c12i = forceFieldParams.iparams[i * (atnr + 1)].lj.c12;
 177             c6j  = forceFieldParams.iparams[j * (atnr + 1)].lj.c6;
 178             c12j = forceFieldParams.iparams[j * (atnr + 1)].lj.c12;
 179             c6   = std::sqrt(c6i * c6j);
 180             if (forceRec.ljpme_combination_rule == LongRangeVdW::LB && !gmx_numzero(c6)
 181                 && !gmx_numzero(c12i) && !gmx_numzero(c12j))
 182             {
 183                 sigmai = gmx::sixthroot(c12i / c6i);
 184                 sigmaj = gmx::sixthroot(c12j / c6j);
 185                 epsi   = c6i * c6i / c12i;
 186                 epsj   = c6j * c6j / c12j;
 187                 c6     = std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
 188             }
 189             /* Store the elements at the same relative positions as C6 in nbfp in order
 190              * to simplify access in the kernels
 191              */
 192             grid[2 * (atnr * i + j)] = c6 * 6.0;
 193         }
 194     }
 195     return grid;
 196 }
 197
 198 enum
 199 {
 200     acNONE = 0,
 201     acCONSTRAINT,
 202     acSETTLE
 203 };
 204
 205 static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t& mtop, const t_forcerec* fr)
 206 {
 207     gmx_bool* type_VDW;
 208     int*      a_con;
 209
 210     snew(type_VDW, fr->ntype);
 211     for (int ai = 0; ai < fr->ntype; ai++)
 212     {
 213         type_VDW[ai] = FALSE;
 214         for (int j = 0; j < fr->ntype; j++)
 215         {
 216             type_VDW[ai] = type_VDW[ai] || fr->bBHAM || C6(fr->nbfp, fr->ntype, ai, j) != 0
 217                            || C12(fr->nbfp, fr->ntype, ai, j) != 0;
 218         }
 219     }
 220
 221     std::vector<cginfo_mb_t> cginfoPerMolblock;
 222     int                      a_offset = 0;
 223     for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
 224     {
 225         const gmx_molblock_t& molb = mtop.molblock[mb];
 226         const gmx_moltype_t&  molt = mtop.moltype[molb.type];
 227         const auto&           excl = molt.excls;
 228
 229         /* Check if the cginfo is identical for all molecules in this block.
 230          * If so, we only need an array of the size of one molecule.
 231          * Otherwise we make an array of #mol times #cgs per molecule.
 232          */
 233         gmx_bool bId = TRUE;
 234         for (int m = 0; m < molb.nmol; m++)
 235         {
 236             const int am = m * molt.atoms.nr;
 237             for (int a = 0; a < molt.atoms.nr; a++)
 238             {
 239                 if (getGroupType(mtop.groups, SimulationAtomGroupType::QuantumMechanics, a_offset + am + a)
 240                     != getGroupType(mtop.groups, SimulationAtomGroupType::QuantumMechanics, a_offset + a))
 241                 {
 242                     bId = FALSE;
 243                 }
 244                 if (!mtop.groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].empty())
 245                 {
 246                     if (mtop.groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics][a_offset + am + a]
 247                         != mtop.groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics][a_offset + a])
 248                     {
 249                         bId = FALSE;
 250                     }
 251                 }
 252             }
 253         }
 254
 255         cginfo_mb_t cginfo_mb;
 256         cginfo_mb.cg_start = a_offset;
 257         cginfo_mb.cg_end   = a_offset + molb.nmol * molt.atoms.nr;
 258         cginfo_mb.cg_mod   = (bId ? 1 : molb.nmol) * molt.atoms.nr;
 259         cginfo_mb.cginfo.resize(cginfo_mb.cg_mod);
 260         gmx::ArrayRef<int> cginfo = cginfo_mb.cginfo;
 261
 262         /* Set constraints flags for constrained atoms */
 263         snew(a_con, molt.atoms.nr);
 264         for (int ftype = 0; ftype < F_NRE; ftype++)
 265         {
 266             if (interaction_function[ftype].flags & IF_CONSTRAINT)
 267             {
 268                 const int nral = NRAL(ftype);
 269                 for (int ia = 0; ia < molt.ilist[ftype].size(); ia += 1 + nral)
 270                 {
 271                     int a;
 272
 273                     for (a = 0; a < nral; a++)
 274                     {
 275                         a_con[molt.ilist[ftype].iatoms[ia + 1 + a]] =
 276                                 (ftype == F_SETTLE ? acSETTLE : acCONSTRAINT);
 277                     }
 278                 }
 279             }
 280         }
 281
 282         for (int m = 0; m < (bId ? 1 : molb.nmol); m++)
 283         {
 284             const int molculeOffsetInBlock = m * molt.atoms.nr;
 285             for (int a = 0; a < molt.atoms.nr; a++)
 286             {
 287                 const t_atom& atom     = molt.atoms.atom[a];
 288                 int&          atomInfo = cginfo[molculeOffsetInBlock + a];
 289
 290                 /* Store the energy group in cginfo */
 291                 int gid = getGroupType(mtop.groups,
 292                                        SimulationAtomGroupType::EnergyOutput,
 293                                        a_offset + molculeOffsetInBlock + a);
 294                 SET_CGINFO_GID(atomInfo, gid);
 295
 296                 bool bHaveVDW = (type_VDW[atom.type] || type_VDW[atom.typeB]);
 297                 bool bHaveQ   = (atom.q != 0 || atom.qB != 0);
 298
 299                 bool haveExclusions = false;
 300                 /* Loop over all the exclusions of atom ai */
 301                 for (const int j : excl[a])
 302                 {
 303                     if (j != a)
 304                     {
 305                         haveExclusions = true;
 306                         break;
 307                     }
 308                 }
 309
 310                 switch (a_con[a])
 311                 {
 312                     case acCONSTRAINT: SET_CGINFO_CONSTR(atomInfo); break;
 313                     case acSETTLE: SET_CGINFO_SETTLE(atomInfo); break;
 314                     default: break;
 315                 }
 316                 if (haveExclusions)
 317                 {
 318                     SET_CGINFO_EXCL_INTER(atomInfo);
 319                 }
 320                 if (bHaveVDW)
 321                 {
 322                     SET_CGINFO_HAS_VDW(atomInfo);
 323                 }
 324                 if (bHaveQ)
 325                 {
 326                     SET_CGINFO_HAS_Q(atomInfo);
 327                 }
 328                 if (fr->efep != FreeEnergyPerturbationType::No && PERTURBED(atom))
 329                 {
 330                     SET_CGINFO_FEP(atomInfo);
 331                 }
 332             }
 333         }
 334
 335         sfree(a_con);
 336
 337         cginfoPerMolblock.push_back(cginfo_mb);
 338
 339         a_offset += molb.nmol * molt.atoms.nr;
 340     }
 341     sfree(type_VDW);
 342
 343     return cginfoPerMolblock;
 344 }
 345
 346 static std::vector<int> cginfo_expand(const int nmb, gmx::ArrayRef<const cginfo_mb_t> cgi_mb)
 347 {
 348     const int ncg = cgi_mb[nmb - 1].cg_end;
 349
 350     std::vector<int> cginfo(ncg);
 351
 352     int mb = 0;
 353     for (int cg = 0; cg < ncg; cg++)
 354     {
 355         while (cg >= cgi_mb[mb].cg_end)
 356         {
 357             mb++;
 358         }
 359         cginfo[cg] = cgi_mb[mb].cginfo[(cg - cgi_mb[mb].cg_start) % cgi_mb[mb].cg_mod];
 360     }
 361
 362     return cginfo;
 363 }
 364
 365 /* Sets the sum of charges (squared) and C6 in the system in fr.
 366  * Returns whether the system has a net charge.
 367  */
 368 static bool set_chargesum(FILE* log, t_forcerec* fr, const gmx_mtop_t& mtop)
 369 {
 370     /*This now calculates sum for q and c6*/
 371     double qsum, q2sum, q, c6sum, c6;
 372
 373     qsum  = 0;
 374     q2sum = 0;
 375     c6sum = 0;
 376     for (const gmx_molblock_t& molb : mtop.molblock)
 377     {
 378         int            nmol  = molb.nmol;
 379         const t_atoms* atoms = &mtop.moltype[molb.type].atoms;
 380         for (int i = 0; i < atoms->nr; i++)
 381         {
 382             q = atoms->atom[i].q;
 383             qsum += nmol * q;
 384             q2sum += nmol * q * q;
 385             c6 = mtop.ffparams.iparams[atoms->atom[i].type * (mtop.ffparams.atnr + 1)].lj.c6;
 386             c6sum += nmol * c6;
 387         }
 388     }
 389     fr->qsum[0]  = qsum;
 390     fr->q2sum[0] = q2sum;
 391     fr->c6sum[0] = c6sum;
 392
 393     if (fr->efep != FreeEnergyPerturbationType::No)
 394     {
 395         qsum  = 0;
 396         q2sum = 0;
 397         c6sum = 0;
 398         for (const gmx_molblock_t& molb : mtop.molblock)
 399         {
 400             int            nmol  = molb.nmol;
 401             const t_atoms* atoms = &mtop.moltype[molb.type].atoms;
 402             for (int i = 0; i < atoms->nr; i++)
 403             {
 404                 q = atoms->atom[i].qB;
 405                 qsum += nmol * q;
 406                 q2sum += nmol * q * q;
 407                 c6 = mtop.ffparams.iparams[atoms->atom[i].typeB * (mtop.ffparams.atnr + 1)].lj.c6;
 408                 c6sum += nmol * c6;
 409             }
 410             fr->qsum[1]  = qsum;
 411             fr->q2sum[1] = q2sum;
 412             fr->c6sum[1] = c6sum;
 413         }
 414     }
 415     else
 416     {
 417         fr->qsum[1]  = fr->qsum[0];
 418         fr->q2sum[1] = fr->q2sum[0];
 419         fr->c6sum[1] = fr->c6sum[0];
 420     }
 421     if (log)
 422     {
 423         if (fr->efep == FreeEnergyPerturbationType::No)
 424         {
 425             fprintf(log, "System total charge: %.3f\n", fr->qsum[0]);
 426         }
 427         else
 428         {
 429             fprintf(log, "System total charge, top. A: %.3f top. B: %.3f\n", fr->qsum[0], fr->qsum[1]);
 430         }
 431     }
 432
 433     /* A cut-off of 1e-4 is used to catch rounding errors due to ascii input */
 434     return (std::abs(fr->qsum[0]) > 1e-4 || std::abs(fr->qsum[1]) > 1e-4);
 435 }
 436
 437 static real calcBuckinghamBMax(FILE* fplog, const gmx_mtop_t& mtop)
 438 {
 439     const t_atoms *at1, *at2;
 440     int            i, j, tpi, tpj, ntypes;
 441     real           b, bmin;
 442
 443     if (fplog)
 444     {
 445         fprintf(fplog, "Determining largest Buckingham b parameter for table\n");
 446     }
 447     ntypes = mtop.ffparams.atnr;
 448
 449     bmin            = -1;
 450     real bham_b_max = 0;
 451     for (size_t mt1 = 0; mt1 < mtop.moltype.size(); mt1++)
 452     {
 453         at1 = &mtop.moltype[mt1].atoms;
 454         for (i = 0; (i < at1->nr); i++)
 455         {
 456             tpi = at1->atom[i].type;
 457             if (tpi >= ntypes)
 458             {
 459                 gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", i, tpi, ntypes);
 460             }
 461
 462             for (size_t mt2 = mt1; mt2 < mtop.moltype.size(); mt2++)
 463             {
 464                 at2 = &mtop.moltype[mt2].atoms;
 465                 for (j = 0; (j < at2->nr); j++)
 466                 {
 467                     tpj = at2->atom[j].type;
 468                     if (tpj >= ntypes)
 469                     {
 470                         gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", j, tpj, ntypes);
 471                     }
 472                     b = mtop.ffparams.iparams[tpi * ntypes + tpj].bham.b;
 473                     if (b > bham_b_max)
 474                     {
 475                         bham_b_max = b;
 476                     }
 477                     if ((b < bmin) || (bmin == -1))
 478                     {
 479                         bmin = b;
 480                     }
 481                 }
 482             }
 483         }
 484     }
 485     if (fplog)
 486     {
 487         fprintf(fplog, "Buckingham b parameters, min: %g, max: %g\n", bmin, bham_b_max);
 488     }
 489
 490     return bham_b_max;
 491 }
 492
 493 /*!\brief If there's bonded interactions of type \c ftype1 or \c
 494  * ftype2 present in the topology, build an array of the number of
 495  * interactions present for each bonded interaction index found in the
 496  * topology.
 497  *
 498  * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
 499  * valid type with that parameter.
 500  *
 501  * \c count will be reallocated as necessary to fit the largest bonded
 502  * interaction index found, and its current size will be returned in
 503  * \c ncount. It will contain zero for every bonded interaction index
 504  * for which no interactions are present in the topology.
 505  */
 506 static void count_tables(int ftype1, int ftype2, const gmx_mtop_t& mtop, int* ncount, int** count)
 507 {
 508     int ftype, i, j, tabnr;
 509
 510     // Loop over all moleculetypes
 511     for (const gmx_moltype_t& molt : mtop.moltype)
 512     {
 513         // Loop over all interaction types
 514         for (ftype = 0; ftype < F_NRE; ftype++)
 515         {
 516             // If the current interaction type is one of the types whose tables we're trying to count...
 517             if (ftype == ftype1 || ftype == ftype2)
 518             {
 519                 const InteractionList& il     = molt.ilist[ftype];
 520                 const int              stride = 1 + NRAL(ftype);
 521                 // ... and there are actually some interactions for this type
 522                 for (i = 0; i < il.size(); i += stride)
 523                 {
 524                     // Find out which table index the user wanted
 525                     tabnr = mtop.ffparams.iparams[il.iatoms[i]].tab.table;
 526                     if (tabnr < 0)
 527                     {
 528                         gmx_fatal(FARGS, "A bonded table number is smaller than 0: %d\n", tabnr);
 529                     }
 530                     // Make room for this index in the data structure
 531                     if (tabnr >= *ncount)
 532                     {
 533                         srenew(*count, tabnr + 1);
 534                         for (j = *ncount; j < tabnr + 1; j++)
 535                         {
 536                             (*count)[j] = 0;
 537                         }
 538                         *ncount = tabnr + 1;
 539                     }
 540                     // Record that this table index is used and must have a valid file
 541                     (*count)[tabnr]++;
 542                 }
 543             }
 544         }
 545     }
 546 }
 547
 548 /*!\brief If there's bonded interactions of flavour \c tabext and type
 549  * \c ftype1 or \c ftype2 present in the topology, seek them in the
 550  * list of filenames passed to mdrun, and make bonded tables from
 551  * those files.
 552  *
 553  * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
 554  * valid type with that parameter.
 555  *
 556  * A fatal error occurs if no matching filename is found.
 557  */
 558 static std::vector<bondedtable_t> make_bonded_tables(FILE*                            fplog,
 559                                                      int                              ftype1,
 560                                                      int                              ftype2,
 561                                                      const gmx_mtop_t&                mtop,
 562                                                      gmx::ArrayRef<const std::string> tabbfnm,
 563                                                      const char*                      tabext)
 564 {
 565     std::vector<bondedtable_t> tab;
 566
 567     int  ncount = 0;
 568     int* count  = nullptr;
 569     count_tables(ftype1, ftype2, mtop, &ncount, &count);
 570
 571     // Are there any relevant tabulated bond interactions?
 572     if (ncount > 0)
 573     {
 574         tab.resize(ncount);
 575         for (int i = 0; i < ncount; i++)
 576         {
 577             // Do any interactions exist that requires this table?
 578             if (count[i] > 0)
 579             {
 580                 // This pattern enforces the current requirement that
 581                 // table filenames end in a characteristic sequence
 582                 // before the file type extension, and avoids table 13
 583                 // being recognized and used for table 1.
 584                 std::string patternToFind = gmx::formatString("_%s%d.%s", tabext, i, ftp2ext(efXVG));
 585                 bool        madeTable     = false;
 586                 for (gmx::index j = 0; j < tabbfnm.ssize() && !madeTable; ++j)
 587                 {
 588                     if (gmx::endsWith(tabbfnm[j], patternToFind))
 589                     {
 590                         // Finally read the table from the file found
 591                         tab[i]    = make_bonded_table(fplog, tabbfnm[j].c_str(), NRAL(ftype1) - 2);
 592                         madeTable = true;
 593                     }
 594                 }
 595                 if (!madeTable)
 596                 {
 597                     bool isPlural = (ftype2 != -1);
 598                     gmx_fatal(FARGS,
 599                               "Tabulated interaction of type '%s%s%s' with index %d cannot be used "
 600                               "because no table file whose name matched '%s' was passed via the "
 601                               "gmx mdrun -tableb command-line option.",
 602                               interaction_function[ftype1].longname,
 603                               isPlural ? "' or '" : "",
 604                               isPlural ? interaction_function[ftype2].longname : "",
 605                               i,
 606                               patternToFind.c_str());
 607                 }
 608             }
 609         }
 610         sfree(count);
 611     }
 612
 613     return tab;
 614 }
 615
 616 void forcerec_set_ranges(t_forcerec* fr, int natoms_force, int natoms_force_constr, int natoms_f_novirsum)
 617 {
 618     fr->natoms_force        = natoms_force;
 619     fr->natoms_force_constr = natoms_force_constr;
 620
 621     for (auto& forceHelperBuffers : fr->forceHelperBuffers)
 622     {
 623         forceHelperBuffers.resize(natoms_f_novirsum);
 624     }
 625 }
 626
 627 static real cutoff_inf(real cutoff)
 628 {
 629     if (cutoff == 0)
 630     {
 631         cutoff = GMX_CUTOFF_INF;
 632     }
 633
 634     return cutoff;
 635 }
 636
 637 /*! \brief Print Coulomb Ewald citations and set ewald coefficients */
 638 static void initCoulombEwaldParameters(FILE*                fp,
 639                                        const t_inputrec&    ir,
 640                                        bool                 systemHasNetCharge,
 641                                        interaction_const_t* ic)
 642 {
 643     if (!EEL_PME_EWALD(ir.coulombtype))
 644     {
 645         return;
 646     }
 647
 648     if (fp)
 649     {
 650         fprintf(fp, "Will do PME sum in reciprocal space for electrostatic interactions.\n");
 651
 652         if (ir.coulombtype == CoulombInteractionType::P3mAD)
 653         {
 654             please_cite(fp, "Hockney1988");
 655             please_cite(fp, "Ballenegger2012");
 656         }
 657         else
 658         {
 659             please_cite(fp, "Essmann95a");
 660         }
 661
 662         if (ir.ewald_geometry == EwaldGeometry::ThreeDC)
 663         {
 664             if (fp)
 665             {
 666                 fprintf(fp,
 667                         "Using the Ewald3DC correction for systems with a slab geometry%s.\n",
 668                         systemHasNetCharge ? " and net charge" : "");
 669             }
 670             please_cite(fp, "In-Chul99a");
 671             if (systemHasNetCharge)
 672             {
 673                 please_cite(fp, "Ballenegger2009");
 674             }
 675         }
 676     }
 677
 678     ic->ewaldcoeff_q = calc_ewaldcoeff_q(ir.rcoulomb, ir.ewald_rtol);
 679     if (fp)
 680     {
 681         fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for Ewald\n", 1 / ic->ewaldcoeff_q);
 682     }
 683
 684     if (ic->coulomb_modifier == InteractionModifiers::PotShift)
 685     {
 686         GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero");
 687         ic->sh_ewald = std::erfc(ic->ewaldcoeff_q * ic->rcoulomb) / ic->rcoulomb;
 688     }
 689     else
 690     {
 691         ic->sh_ewald = 0;
 692     }
 693 }
 694
 695 /*! \brief Print Van der Waals Ewald citations and set ewald coefficients */
 696 static void initVdwEwaldParameters(FILE* fp, const t_inputrec& ir, interaction_const_t* ic)
 697 {
 698     if (!EVDW_PME(ir.vdwtype))
 699     {
 700         return;
 701     }
 702
 703     if (fp)
 704     {
 705         fprintf(fp, "Will do PME sum in reciprocal space for LJ dispersion interactions.\n");
 706         please_cite(fp, "Essmann95a");
 707     }
 708     ic->ewaldcoeff_lj = calc_ewaldcoeff_lj(ir.rvdw, ir.ewald_rtol_lj);
 709     if (fp)
 710     {
 711         fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for LJ Ewald\n", 1 / ic->ewaldcoeff_lj);
 712     }
 713
 714     if (ic->vdw_modifier == InteractionModifiers::PotShift)
 715     {
 716         real crc2       = gmx::square(ic->ewaldcoeff_lj * ic->rvdw);
 717         ic->sh_lj_ewald = (std::exp(-crc2) * (1 + crc2 + 0.5 * crc2 * crc2) - 1) / gmx::power6(ic->rvdw);
 718     }
 719     else
 720     {
 721         ic->sh_lj_ewald = 0;
 722     }
 723 }
 724
 725 /* Generate Coulomb and/or Van der Waals Ewald long-range correction tables
 726  *
 727  * Tables are generated for one or both, depending on if the pointers
 728  * are non-null. The spacing for both table sets is the same and obeys
 729  * both accuracy requirements, when relevant.
 730  */
 731 static void init_ewald_f_table(const interaction_const_t& ic,
 732                                const real                 rlist,
 733                                const real                 tabext,
 734                                EwaldCorrectionTables*     coulombTables,
 735                                EwaldCorrectionTables*     vdwTables)
 736 {
 737     const bool useCoulombTable = (EEL_PME_EWALD(ic.eeltype) && coulombTables != nullptr);
 738     const bool useVdwTable     = (EVDW_PME(ic.vdwtype) && vdwTables != nullptr);
 739
 740     /* Get the Ewald table spacing based on Coulomb and/or LJ
 741      * Ewald coefficients and rtol.
 742      */
 743     const real tableScale = ewald_spline3_table_scale(ic, useCoulombTable, useVdwTable);
 744
 745     const bool havePerturbedNonbondeds = (ic.softCoreParameters != nullptr);
 746
 747     real tableLen = ic.rcoulomb;
 748     if ((useCoulombTable || useVdwTable) && havePerturbedNonbondeds && rlist + tabext > 0.0)
 749     {
 750         /* TODO: Ideally this should also check if couple-intramol == no, but that isn't
 751          * stored in ir. Grompp puts that info into an opts structure that doesn't make it into the tpr.
 752          * The alternative is to look through all the exclusions and check if they come from
 753          * couple-intramol == no. Meanwhile, always having larger tables should only affect
 754          * memory consumption, not speed (barring cache issues).
 755          */
 756         tableLen = rlist + tabext;
 757     }
 758     const int tableSize = static_cast<int>(tableLen * tableScale) + 2;
 759
 760     if (useCoulombTable)
 761     {
 762         *coulombTables =
 763                 generateEwaldCorrectionTables(tableSize, tableScale, ic.ewaldcoeff_q, v_q_ewald_lr);
 764     }
 765
 766     if (useVdwTable)
 767     {
 768         *vdwTables = generateEwaldCorrectionTables(tableSize, tableScale, ic.ewaldcoeff_lj, v_lj_ewald_lr);
 769     }
 770 }
 771
 772 void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real rlist, const real tableExtensionLength)
 773 {
 774     if (EEL_PME_EWALD(ic->eeltype) || EVDW_PME(ic->vdwtype))
 775     {
 776         init_ewald_f_table(
 777                 *ic, rlist, tableExtensionLength, ic->coulombEwaldTables.get(), ic->vdwEwaldTables.get());
 778         if (fp != nullptr)
 779         {
 780             fprintf(fp,
 781                     "Initialized non-bonded Ewald tables, spacing: %.2e size: %zu\n\n",
 782                     1 / ic->coulombEwaldTables->scale,
 783                     ic->coulombEwaldTables->tableF.size());
 784         }
 785     }
 786 }
 787
 788 static void clear_force_switch_constants(shift_consts_t* sc)
 789 {
 790     sc->c2   = 0;
 791     sc->c3   = 0;
 792     sc->cpot = 0;
 793 }
 794
 795 static void force_switch_constants(real p, real rsw, real rc, shift_consts_t* sc)
 796 {
 797     /* Here we determine the coefficient for shifting the force to zero
 798      * between distance rsw and the cut-off rc.
 799      * For a potential of r^-p, we have force p*r^-(p+1).
 800      * But to save flops we absorb p in the coefficient.
 801      * Thus we get:
 802      * force/p   = r^-(p+1) + c2*r^2 + c3*r^3
 803      * potential = r^-p + c2/3*r^3 + c3/4*r^4 + cpot
 804      */
 805     sc->c2   = ((p + 1) * rsw - (p + 4) * rc) / (pow(rc, p + 2) * gmx::square(rc - rsw));
 806     sc->c3   = -((p + 1) * rsw - (p + 3) * rc) / (pow(rc, p + 2) * gmx::power3(rc - rsw));
 807     sc->cpot = -pow(rc, -p) + p * sc->c2 / 3 * gmx::power3(rc - rsw)
 808                + p * sc->c3 / 4 * gmx::power4(rc - rsw);
 809 }
 810
 811 static void potential_switch_constants(real rsw, real rc, switch_consts_t* sc)
 812 {
 813     /* The switch function is 1 at rsw and 0 at rc.
 814      * The derivative and second derivate are zero at both ends.
 815      * rsw        = max(r - r_switch, 0)
 816      * sw         = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5
 817      * dsw        = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4
 818      * force      = force*dsw - potential*sw
 819      * potential *= sw
 820      */
 821     sc->c3 = -10 / gmx::power3(rc - rsw);
 822     sc->c4 = 15 / gmx::power4(rc - rsw);
 823     sc->c5 = -6 / gmx::power5(rc - rsw);
 824 }
 825
 826 /*! \brief Construct interaction constants
 827  *
 828  * This data is used (particularly) by search and force code for
 829  * short-range interactions. Many of these are constant for the whole
 830  * simulation; some are constant only after PME tuning completes.
 831  */
 832 static void init_interaction_const(FILE*                 fp,
 833                                    interaction_const_t** interaction_const,
 834                                    const t_inputrec&     ir,
 835                                    const gmx_mtop_t&     mtop,
 836                                    bool                  systemHasNetCharge)
 837 {
 838     interaction_const_t* ic = new interaction_const_t;
 839
 840     ic->coulombEwaldTables = std::make_unique<EwaldCorrectionTables>();
 841     ic->vdwEwaldTables     = std::make_unique<EwaldCorrectionTables>();
 842
 843     /* Lennard-Jones */
 844     ic->vdwtype         = ir.vdwtype;
 845     ic->vdw_modifier    = ir.vdw_modifier;
 846     ic->reppow          = mtop.ffparams.reppow;
 847     ic->rvdw            = cutoff_inf(ir.rvdw);
 848     ic->rvdw_switch     = ir.rvdw_switch;
 849     ic->ljpme_comb_rule = ir.ljpme_combination_rule;
 850     ic->useBuckingham   = (mtop.ffparams.functype[0] == F_BHAM);
 851     if (ic->useBuckingham)
 852     {
 853         ic->buckinghamBMax = calcBuckinghamBMax(fp, mtop);
 854     }
 855
 856     initVdwEwaldParameters(fp, ir, ic);
 857
 858     clear_force_switch_constants(&ic->dispersion_shift);
 859     clear_force_switch_constants(&ic->repulsion_shift);
 860
 861     switch (ic->vdw_modifier)
 862     {
 863         case InteractionModifiers::PotShift:
 864             /* Only shift the potential, don't touch the force */
 865             ic->dispersion_shift.cpot = -1.0 / gmx::power6(ic->rvdw);
 866             ic->repulsion_shift.cpot  = -1.0 / gmx::power12(ic->rvdw);
 867             break;
 868         case InteractionModifiers::ForceSwitch:
 869             /* Switch the force, switch and shift the potential */
 870             force_switch_constants(6.0, ic->rvdw_switch, ic->rvdw, &ic->dispersion_shift);
 871             force_switch_constants(12.0, ic->rvdw_switch, ic->rvdw, &ic->repulsion_shift);
 872             break;
 873         case InteractionModifiers::PotSwitch:
 874             /* Switch the potential and force */
 875             potential_switch_constants(ic->rvdw_switch, ic->rvdw, &ic->vdw_switch);
 876             break;
 877         case InteractionModifiers::None:
 878         case InteractionModifiers::ExactCutoff:
 879             /* Nothing to do here */
 880             break;
 881         default: gmx_incons("unimplemented potential modifier");
 882     }
 883
 884     /* Electrostatics */
 885     ic->eeltype          = ir.coulombtype;
 886     ic->coulomb_modifier = ir.coulomb_modifier;
 887     ic->rcoulomb         = cutoff_inf(ir.rcoulomb);
 888     ic->rcoulomb_switch  = ir.rcoulomb_switch;
 889     ic->epsilon_r        = ir.epsilon_r;
 890
 891     /* Set the Coulomb energy conversion factor */
 892     if (ic->epsilon_r != 0)
 893     {
 894         ic->epsfac = ONE_4PI_EPS0 / ic->epsilon_r;
 895     }
 896     else
 897     {
 898         /* eps = 0 is infinite dieletric: no Coulomb interactions */
 899         ic->epsfac = 0;
 900     }
 901
 902     /* Reaction-field */
 903     if (EEL_RF(ic->eeltype))
 904     {
 905         GMX_RELEASE_ASSERT(ic->eeltype != CoulombInteractionType::GRFNotused,
 906                            "GRF is no longer supported");
 907         ic->reactionFieldPermitivity = ir.epsilon_rf;
 908
 909         calc_rffac(fp,
 910                    ic->epsilon_r,
 911                    ic->reactionFieldPermitivity,
 912                    ic->rcoulomb,
 913                    &ic->reactionFieldCoefficient,
 914                    &ic->reactionFieldShift);
 915     }
 916     else
 917     {
 918         /* For plain cut-off we might use the reaction-field kernels */
 919         ic->reactionFieldPermitivity = ic->epsilon_r;
 920         ic->reactionFieldCoefficient = 0;
 921         if (ir.coulomb_modifier == InteractionModifiers::PotShift)
 922         {
 923             ic->reactionFieldShift = 1 / ic->rcoulomb;
 924         }
 925         else
 926         {
 927             ic->reactionFieldShift = 0;
 928         }
 929     }
 930
 931     initCoulombEwaldParameters(fp, ir, systemHasNetCharge, ic);
 932
 933     if (fp != nullptr)
 934     {
 935         real dispersion_shift;
 936
 937         dispersion_shift = ic->dispersion_shift.cpot;
 938         if (EVDW_PME(ic->vdwtype))
 939         {
 940             dispersion_shift -= ic->sh_lj_ewald;
 941         }
 942         fprintf(fp, "Potential shift: LJ r^-12: %.3e r^-6: %.3e", ic->repulsion_shift.cpot, dispersion_shift);
 943
 944         if (ic->eeltype == CoulombInteractionType::Cut)
 945         {
 946             fprintf(fp, ", Coulomb %.e", -ic->reactionFieldShift);
 947         }
 948         else if (EEL_PME(ic->eeltype))
 949         {
 950             fprintf(fp, ", Ewald %.3e", -ic->sh_ewald);
 951         }
 952         fprintf(fp, "\n");
 953     }
 954
 955     if (ir.efep != FreeEnergyPerturbationType::No)
 956     {
 957         GMX_RELEASE_ASSERT(ir.fepvals, "ir.fepvals should be set wth free-energy");
 958         ic->softCoreParameters = std::make_unique<interaction_const_t::SoftCoreParameters>(*ir.fepvals);
 959     }
 960
 961     *interaction_const = ic;
 962 }
 963
 964 void init_forcerec(FILE*                            fp,
 965                    const gmx::MDLogger&             mdlog,
 966                    t_forcerec*                      fr,
 967                    const t_inputrec&                ir,
 968                    const gmx_mtop_t&                mtop,
 969                    const t_commrec*                 cr,
 970                    matrix                           box,
 971                    const char*                      tabfn,
 972                    const char*                      tabpfn,
 973                    gmx::ArrayRef<const std::string> tabbfnm,
 974                    real                             print_force)
 975 {
 976     /* The CMake default turns SIMD kernels on, but it might be turned off further down... */
 977     fr->use_simd_kernels = GMX_USE_SIMD_KERNELS;
 978
 979     if (check_box(ir.pbcType, box))
 980     {
 981         gmx_fatal(FARGS, "%s", check_box(ir.pbcType, box));
 982     }
 983
 984     /* Test particle insertion ? */
 985     if (EI_TPI(ir.eI))
 986     {
 987         /* Set to the size of the molecule to be inserted (the last one) */
 988         gmx::RangePartitioning molecules = gmx_mtop_molecules(mtop);
 989         fr->n_tpi                        = molecules.block(molecules.numBlocks() - 1).size();
 990     }
 991     else
 992     {
 993         fr->n_tpi = 0;
 994     }
 995
 996     if (ir.coulombtype == CoulombInteractionType::RFNecUnsupported
 997         || ir.coulombtype == CoulombInteractionType::GRFNotused)
 998     {
 999         gmx_fatal(FARGS, "%s electrostatics is no longer supported", enumValueToString(ir.coulombtype));
1000     }
1001
1002     if (ir.bAdress)
1003     {
1004         gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1005     }
1006     if (ir.useTwinRange)
1007     {
1008         gmx_fatal(FARGS, "Twin-range simulations are no longer supported");
1009     }
1010     /* Copy the user determined parameters */
1011     fr->userint1  = ir.userint1;
1012     fr->userint2  = ir.userint2;
1013     fr->userint3  = ir.userint3;
1014     fr->userint4  = ir.userint4;
1015     fr->userreal1 = ir.userreal1;
1016     fr->userreal2 = ir.userreal2;
1017     fr->userreal3 = ir.userreal3;
1018     fr->userreal4 = ir.userreal4;
1019
1020     /* Shell stuff */
1021     fr->fc_stepsize = ir.fc_stepsize;
1022
1023     /* Free energy */
1024     fr->efep = ir.efep;
1025
1026     if ((getenv("GMX_DISABLE_SIMD_KERNELS") != nullptr) || (getenv("GMX_NOOPTIMIZEDKERNELS") != nullptr))
1027     {
1028         fr->use_simd_kernels = FALSE;
1029         if (fp != nullptr)
1030         {
1031             fprintf(fp,
1032                     "\nFound environment variable GMX_DISABLE_SIMD_KERNELS.\n"
1033                     "Disabling the usage of any SIMD-specific non-bonded & bonded kernel routines\n"
1034                     "(e.g. SSE2/SSE4.1/AVX).\n\n");
1035         }
1036     }
1037
1038     fr->bBHAM = (mtop.ffparams.functype[0] == F_BHAM);
1039
1040     /* Neighbour searching stuff */
1041     fr->pbcType = ir.pbcType;
1042
1043     /* Determine if we will do PBC for distances in bonded interactions */
1044     if (fr->pbcType == PbcType::No)
1045     {
1046         fr->bMolPBC = FALSE;
1047     }
1048     else
1049     {
1050         const bool useEwaldSurfaceCorrection = (EEL_PME_EWALD(ir.coulombtype) && ir.epsilon_surface != 0);
1051         const bool haveOrientationRestraints = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
1052         if (!DOMAINDECOMP(cr))
1053         {
1054             fr->bMolPBC = true;
1055
1056             if (useEwaldSurfaceCorrection || haveOrientationRestraints)
1057             {
1058                 fr->wholeMoleculeTransform =
1059                         std::make_unique<gmx::WholeMoleculeTransform>(mtop, ir.pbcType);
1060             }
1061         }
1062         else
1063         {
1064             fr->bMolPBC = dd_bonded_molpbc(cr->dd, fr->pbcType);
1065
1066             /* With Ewald surface correction it is useful to support e.g. running water
1067              * in parallel with update groups.
1068              * With orientation restraints there is no sensible use case supported with DD.
1069              */
1070             if ((useEwaldSurfaceCorrection && !dd_moleculesAreAlwaysWhole(*cr->dd)) || haveOrientationRestraints)
1071             {
1072                 gmx_fatal(FARGS,
1073                           "You requested Ewald surface correction or orientation restraints, "
1074                           "but molecules are broken "
1075                           "over periodic boundary conditions by the domain decomposition. "
1076                           "Run without domain decomposition instead.");
1077             }
1078         }
1079
1080         if (useEwaldSurfaceCorrection)
1081         {
1082             GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr) || dd_moleculesAreAlwaysWhole(*cr->dd),
1083                                "Molecules can not be broken by PBC with epsilon_surface > 0");
1084         }
1085     }
1086
1087     fr->rc_scaling = ir.refcoord_scaling;
1088     copy_rvec(ir.posres_com, fr->posres_com);
1089     copy_rvec(ir.posres_comB, fr->posres_comB);
1090     fr->rlist                  = cutoff_inf(ir.rlist);
1091     fr->ljpme_combination_rule = ir.ljpme_combination_rule;
1092
1093     /* This now calculates sum for q and c6*/
1094     bool systemHasNetCharge = set_chargesum(fp, fr, mtop);
1095
1096     /* fr->ic is used both by verlet and group kernels (to some extent) now */
1097     init_interaction_const(fp, &fr->ic, ir, mtop, systemHasNetCharge);
1098     init_interaction_const_tables(fp, fr->ic, fr->rlist, ir.tabext);
1099
1100     const interaction_const_t* ic = fr->ic;
1101
1102     /* TODO: Replace this Ewald table or move it into interaction_const_t */
1103     if (ir.coulombtype == CoulombInteractionType::Ewald)
1104     {
1105         init_ewald_tab(&(fr->ewald_table), ir, fp);
1106     }
1107
1108     /* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */
1109     switch (ic->eeltype)
1110     {
1111         case CoulombInteractionType::Cut:
1112             fr->nbkernel_elec_interaction = NbkernelElecType::Coulomb;
1113             break;
1114
1115         case CoulombInteractionType::RF:
1116         case CoulombInteractionType::RFZero:
1117             fr->nbkernel_elec_interaction = NbkernelElecType::ReactionField;
1118             break;
1119
1120         case CoulombInteractionType::Switch:
1121         case CoulombInteractionType::Shift:
1122         case CoulombInteractionType::User:
1123         case CoulombInteractionType::PmeSwitch:
1124         case CoulombInteractionType::PmeUser:
1125         case CoulombInteractionType::PmeUserSwitch:
1126             fr->nbkernel_elec_interaction = NbkernelElecType::CubicSplineTable;
1127             break;
1128
1129         case CoulombInteractionType::Pme:
1130         case CoulombInteractionType::P3mAD:
1131         case CoulombInteractionType::Ewald:
1132             fr->nbkernel_elec_interaction = NbkernelElecType::Ewald;
1133             break;
1134
1135         default:
1136             gmx_fatal(FARGS, "Unsupported electrostatic interaction: %s", enumValueToString(ic->eeltype));
1137     }
1138     fr->nbkernel_elec_modifier = ic->coulomb_modifier;
1139
1140     /* Vdw: Translate from mdp settings to kernel format */
1141     switch (ic->vdwtype)
1142     {
1143         case VanDerWaalsType::Cut:
1144             if (fr->bBHAM)
1145             {
1146                 fr->nbkernel_vdw_interaction = NbkernelVdwType::Buckingham;
1147             }
1148             else
1149             {
1150                 fr->nbkernel_vdw_interaction = NbkernelVdwType::LennardJones;
1151             }
1152             break;
1153         case VanDerWaalsType::Pme: fr->nbkernel_vdw_interaction = NbkernelVdwType::LJEwald; break;
1154
1155         case VanDerWaalsType::Switch:
1156         case VanDerWaalsType::Shift:
1157         case VanDerWaalsType::User:
1158             fr->nbkernel_vdw_interaction = NbkernelVdwType::CubicSplineTable;
1159             break;
1160
1161         default:
1162             gmx_fatal(FARGS, "Unsupported vdw interaction: %s", enumValueToString(ic->vdwtype));
1163     }
1164     fr->nbkernel_vdw_modifier = ic->vdw_modifier;
1165
1166     if (!gmx_within_tol(ic->reppow, 12.0, 10 * GMX_DOUBLE_EPS))
1167     {
1168         gmx_fatal(FARGS, "Only LJ repulsion power 12 is supported");
1169     }
1170     /* Older tpr files can contain Coulomb user tables with the Verlet cutoff-scheme,
1171      * while mdrun does not (and never did) support this.
1172      */
1173     if (EEL_USER(fr->ic->eeltype))
1174     {
1175         gmx_fatal(FARGS, "Electrostatics type %s is currently not supported", enumValueToString(ir.coulombtype));
1176     }
1177
1178     fr->bvdwtab  = FALSE;
1179     fr->bcoultab = FALSE;
1180
1181     /* 1-4 interaction electrostatics */
1182     fr->fudgeQQ = mtop.ffparams.fudgeQQ;
1183
1184     // Multiple time stepping
1185     fr->useMts = ir.useMts;
1186
1187     if (fr->useMts)
1188     {
1189         GMX_ASSERT(gmx::checkMtsRequirements(ir).empty(), "All MTS requirements should be met here");
1190     }
1191
1192     const bool haveDirectVirialContributionsFast = fr->forceProviders->hasForceProvider()
1193                                                    || gmx_mtop_ftype_count(mtop, F_POSRES) > 0
1194                                                    || gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0
1195                                                    || ir.nwall > 0 || ir.bPull || ir.bRot || ir.bIMD;
1196     const bool haveDirectVirialContributionsSlow = EEL_FULL(ic->eeltype) || EVDW_PME(ic->vdwtype);
1197     for (int i = 0; i < (fr->useMts ? 2 : 1); i++)
1198     {
1199         bool haveDirectVirialContributions =
1200                 (((!fr->useMts || i == 0) && haveDirectVirialContributionsFast)
1201                  || ((!fr->useMts || i == 1) && haveDirectVirialContributionsSlow));
1202         fr->forceHelperBuffers.emplace_back(haveDirectVirialContributions);
1203     }
1204
1205     if (fr->shift_vec == nullptr)
1206     {
1207         snew(fr->shift_vec, SHIFTS);
1208     }
1209
1210     if (fr->nbfp.empty())
1211     {
1212         fr->ntype = mtop.ffparams.atnr;
1213         fr->nbfp  = makeNonBondedParameterLists(mtop.ffparams, fr->bBHAM);
1214         if (EVDW_PME(ic->vdwtype))
1215         {
1216             fr->ljpme_c6grid = makeLJPmeC6GridCorrectionParameters(mtop.ffparams, *fr);
1217         }
1218     }
1219
1220     /* Copy the energy group exclusions */
1221     fr->egp_flags = ir.opts.egp_flags;
1222
1223     /* Van der Waals stuff */
1224     if ((ic->vdwtype != VanDerWaalsType::Cut) && (ic->vdwtype != VanDerWaalsType::User) && !fr->bBHAM)
1225     {
1226         if (ic->rvdw_switch >= ic->rvdw)
1227         {
1228             gmx_fatal(FARGS, "rvdw_switch (%f) must be < rvdw (%f)", ic->rvdw_switch, ic->rvdw);
1229         }
1230         if (fp)
1231         {
1232             fprintf(fp,
1233                     "Using %s Lennard-Jones, switch between %g and %g nm\n",
1234                     (ic->eeltype == CoulombInteractionType::Switch) ? "switched" : "shifted",
1235                     ic->rvdw_switch,
1236                     ic->rvdw);
1237         }
1238     }
1239
1240     if (fr->bBHAM && EVDW_PME(ic->vdwtype))
1241     {
1242         gmx_fatal(FARGS, "LJ PME not supported with Buckingham");
1243     }
1244
1245     if (fr->bBHAM && (ic->vdwtype == VanDerWaalsType::Shift || ic->vdwtype == VanDerWaalsType::Switch))
1246     {
1247         gmx_fatal(FARGS, "Switch/shift interaction not supported with Buckingham");
1248     }
1249
1250     if (fr->bBHAM)
1251     {
1252         gmx_fatal(FARGS, "The Verlet cutoff-scheme does not (yet) support Buckingham");
1253     }
1254
1255     if (ir.implicit_solvent)
1256     {
1257         gmx_fatal(FARGS, "Implict solvation is no longer supported.");
1258     }
1259
1260
1261     /* This code automatically gives table length tabext without cut-off's,
1262      * in that case grompp should already have checked that we do not need
1263      * normal tables and we only generate tables for 1-4 interactions.
1264      */
1265     real rtab = ir.rlist + ir.tabext;
1266
1267     /* We want to use unmodified tables for 1-4 coulombic
1268      * interactions, so we must in general have an extra set of
1269      * tables. */
1270     if (gmx_mtop_ftype_count(mtop, F_LJ14) > 0 || gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0
1271         || gmx_mtop_ftype_count(mtop, F_LJC_PAIRS_NB) > 0)
1272     {
1273         fr->pairsTable = make_tables(fp, ic, tabpfn, rtab, GMX_MAKETABLES_14ONLY);
1274     }
1275
1276     /* Wall stuff */
1277     fr->nwall = ir.nwall;
1278     if (ir.nwall && ir.wall_type == WallType::Table)
1279     {
1280         make_wall_tables(fp, ir, tabfn, &mtop.groups, fr);
1281     }
1282
1283     fr->fcdata = std::make_unique<t_fcdata>();
1284
1285     if (!tabbfnm.empty())
1286     {
1287         t_fcdata& fcdata = *fr->fcdata;
1288         // Need to catch std::bad_alloc
1289         // TODO Don't need to catch this here, when merging with master branch
1290         try
1291         {
1292             // TODO move these tables into a separate struct and store reference in ListedForces
1293             fcdata.bondtab  = make_bonded_tables(fp, F_TABBONDS, F_TABBONDSNC, mtop, tabbfnm, "b");
1294             fcdata.angletab = make_bonded_tables(fp, F_TABANGLES, -1, mtop, tabbfnm, "a");
1295             fcdata.dihtab   = make_bonded_tables(fp, F_TABDIHS, -1, mtop, tabbfnm, "d");
1296         }
1297         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1298     }
1299     else
1300     {
1301         if (debug)
1302         {
1303             fprintf(debug,
1304                     "No fcdata or table file name passed, can not read table, can not do bonded "
1305                     "interactions\n");
1306         }
1307     }
1308
1309     /* Initialize the thread working data for bonded interactions */
1310     if (fr->useMts)
1311     {
1312         // Add one ListedForces object for each MTS level
1313         bool isFirstLevel = true;
1314         for (const auto& mtsLevel : ir.mtsLevels)
1315         {
1316             ListedForces::InteractionSelection interactionSelection;
1317             const auto&                        forceGroups = mtsLevel.forceGroups;
1318             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Pair)])
1319             {
1320                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Pairs));
1321             }
1322             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Dihedral)])
1323             {
1324                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Dihedrals));
1325             }
1326             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Angle)])
1327             {
1328                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Angles));
1329             }
1330             if (isFirstLevel)
1331             {
1332                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Rest));
1333                 isFirstLevel = false;
1334             }
1335             fr->listedForces.emplace_back(mtop.ffparams,
1336                                           mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1337                                           gmx_omp_nthreads_get(emntBonded),
1338                                           interactionSelection,
1339                                           fp);
1340         }
1341     }
1342     else
1343     {
1344         // Add one ListedForces object with all listed interactions
1345         fr->listedForces.emplace_back(mtop.ffparams,
1346                                       mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1347                                       gmx_omp_nthreads_get(emntBonded),
1348                                       ListedForces::interactionSelectionAll(),
1349                                       fp);
1350     }
1351
1352     // QM/MM initialization if requested
1353     if (ir.bQMMM)
1354     {
1355         gmx_incons("QM/MM was requested, but is no longer available in GROMACS");
1356     }
1357
1358     /* Set all the static charge group info */
1359     fr->cginfo_mb = init_cginfo_mb(mtop, fr);
1360     if (!DOMAINDECOMP(cr))
1361     {
1362         fr->cginfo = cginfo_expand(mtop.molblock.size(), fr->cginfo_mb);
1363     }
1364
1365     if (!DOMAINDECOMP(cr))
1366     {
1367         forcerec_set_ranges(fr, mtop.natoms, mtop.natoms, mtop.natoms);
1368     }
1369
1370     fr->print_force = print_force;
1371
1372     fr->nthread_ewc = gmx_omp_nthreads_get(emntBonded);
1373     snew(fr->ewc_t, fr->nthread_ewc);
1374
1375     if (ir.eDispCorr != DispersionCorrectionType::No)
1376     {
1377         fr->dispersionCorrection = std::make_unique<DispersionCorrection>(
1378                 mtop, ir, fr->bBHAM, fr->ntype, fr->nbfp, *fr->ic, tabfn);
1379         fr->dispersionCorrection->print(mdlog);
1380     }
1381
1382     if (fp != nullptr)
1383     {
1384         /* Here we switch from using mdlog, which prints the newline before
1385          * the paragraph, to our old fprintf logging, which prints the newline
1386          * after the paragraph, so we should add a newline here.
1387          */
1388         fprintf(fp, "\n");
1389     }
1390 }
1391
1392 t_forcerec::t_forcerec() = default;
1393
1394 t_forcerec::~t_forcerec()
1395 {
1396     /* Note: This code will disappear when types are converted to C++ */
1397     sfree(shift_vec);
1398     sfree(ewc_t);
1399 }