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->haveBuckingham || 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 interaction_const_t init_interaction_const(FILE*             fp,
 833                                                   const t_inputrec& ir,
 834                                                   const gmx_mtop_t& mtop,
 835                                                   bool              systemHasNetCharge)
 836 {
 837     interaction_const_t interactionConst;
 838
 839     interactionConst.coulombEwaldTables = std::make_unique<EwaldCorrectionTables>();
 840     interactionConst.vdwEwaldTables     = std::make_unique<EwaldCorrectionTables>();
 841
 842     /* Lennard-Jones */
 843     interactionConst.vdwtype         = ir.vdwtype;
 844     interactionConst.vdw_modifier    = ir.vdw_modifier;
 845     interactionConst.reppow          = mtop.ffparams.reppow;
 846     interactionConst.rvdw            = cutoff_inf(ir.rvdw);
 847     interactionConst.rvdw_switch     = ir.rvdw_switch;
 848     interactionConst.ljpme_comb_rule = ir.ljpme_combination_rule;
 849     interactionConst.useBuckingham   = (mtop.ffparams.functype[0] == F_BHAM);
 850     if (interactionConst.useBuckingham)
 851     {
 852         interactionConst.buckinghamBMax = calcBuckinghamBMax(fp, mtop);
 853     }
 854
 855     initVdwEwaldParameters(fp, ir, &interactionConst);
 856
 857     clear_force_switch_constants(&interactionConst.dispersion_shift);
 858     clear_force_switch_constants(&interactionConst.repulsion_shift);
 859
 860     switch (interactionConst.vdw_modifier)
 861     {
 862         case InteractionModifiers::PotShift:
 863             /* Only shift the potential, don't touch the force */
 864             interactionConst.dispersion_shift.cpot = -1.0 / gmx::power6(interactionConst.rvdw);
 865             interactionConst.repulsion_shift.cpot  = -1.0 / gmx::power12(interactionConst.rvdw);
 866             break;
 867         case InteractionModifiers::ForceSwitch:
 868             /* Switch the force, switch and shift the potential */
 869             force_switch_constants(
 870                     6.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.dispersion_shift);
 871             force_switch_constants(
 872                     12.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.repulsion_shift);
 873             break;
 874         case InteractionModifiers::PotSwitch:
 875             /* Switch the potential and force */
 876             potential_switch_constants(
 877                     interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.vdw_switch);
 878             break;
 879         case InteractionModifiers::None:
 880         case InteractionModifiers::ExactCutoff:
 881             /* Nothing to do here */
 882             break;
 883         default: gmx_incons("unimplemented potential modifier");
 884     }
 885
 886     /* Electrostatics */
 887     interactionConst.eeltype          = ir.coulombtype;
 888     interactionConst.coulomb_modifier = ir.coulomb_modifier;
 889     interactionConst.rcoulomb         = cutoff_inf(ir.rcoulomb);
 890     interactionConst.rcoulomb_switch  = ir.rcoulomb_switch;
 891     interactionConst.epsilon_r        = ir.epsilon_r;
 892
 893     /* Set the Coulomb energy conversion factor */
 894     if (interactionConst.epsilon_r != 0)
 895     {
 896         interactionConst.epsfac = gmx::c_one4PiEps0 / interactionConst.epsilon_r;
 897     }
 898     else
 899     {
 900         /* eps = 0 is infinite dieletric: no Coulomb interactions */
 901         interactionConst.epsfac = 0;
 902     }
 903
 904     /* Reaction-field */
 905     if (EEL_RF(interactionConst.eeltype))
 906     {
 907         GMX_RELEASE_ASSERT(interactionConst.eeltype != CoulombInteractionType::GRFNotused,
 908                            "GRF is no longer supported");
 909         interactionConst.reactionFieldPermitivity = ir.epsilon_rf;
 910         calc_rffac(fp,
 911                    interactionConst.epsilon_r,
 912                    interactionConst.reactionFieldPermitivity,
 913                    interactionConst.rcoulomb,
 914                    &interactionConst.reactionFieldCoefficient,
 915                    &interactionConst.reactionFieldShift);
 916     }
 917     else
 918     {
 919         /* For plain cut-off we might use the reaction-field kernels */
 920         interactionConst.reactionFieldPermitivity = interactionConst.epsilon_r;
 921         interactionConst.reactionFieldCoefficient = 0;
 922         if (ir.coulomb_modifier == InteractionModifiers::PotShift)
 923         {
 924             interactionConst.reactionFieldShift = 1 / interactionConst.rcoulomb;
 925         }
 926         else
 927         {
 928             interactionConst.reactionFieldShift = 0;
 929         }
 930     }
 931
 932     initCoulombEwaldParameters(fp, ir, systemHasNetCharge, &interactionConst);
 933
 934     if (fp != nullptr)
 935     {
 936         real dispersion_shift;
 937
 938         dispersion_shift = interactionConst.dispersion_shift.cpot;
 939         if (EVDW_PME(interactionConst.vdwtype))
 940         {
 941             dispersion_shift -= interactionConst.sh_lj_ewald;
 942         }
 943         fprintf(fp,
 944                 "Potential shift: LJ r^-12: %.3e r^-6: %.3e",
 945                 interactionConst.repulsion_shift.cpot,
 946                 dispersion_shift);
 947
 948         if (interactionConst.eeltype == CoulombInteractionType::Cut)
 949         {
 950             fprintf(fp, ", Coulomb %.e", -interactionConst.reactionFieldShift);
 951         }
 952         else if (EEL_PME(interactionConst.eeltype))
 953         {
 954             fprintf(fp, ", Ewald %.3e", -interactionConst.sh_ewald);
 955         }
 956         fprintf(fp, "\n");
 957     }
 958
 959     if (ir.efep != FreeEnergyPerturbationType::No)
 960     {
 961         GMX_RELEASE_ASSERT(ir.fepvals, "ir.fepvals should be set wth free-energy");
 962         interactionConst.softCoreParameters =
 963                 std::make_unique<interaction_const_t::SoftCoreParameters>(*ir.fepvals);
 964     }
 965
 966     return interactionConst;
 967 }
 968
 969 void init_forcerec(FILE*                            fplog,
 970                    const gmx::MDLogger&             mdlog,
 971                    t_forcerec*                      forcerec,
 972                    const t_inputrec&                inputrec,
 973                    const gmx_mtop_t&                mtop,
 974                    const t_commrec*                 commrec,
 975                    matrix                           box,
 976                    const char*                      tabfn,
 977                    const char*                      tabpfn,
 978                    gmx::ArrayRef<const std::string> tabbfnm,
 979                    real                             print_force)
 980 {
 981     /* The CMake default turns SIMD kernels on, but it might be turned off further down... */
 982     forcerec->use_simd_kernels = GMX_USE_SIMD_KERNELS;
 983
 984     if (check_box(inputrec.pbcType, box))
 985     {
 986         gmx_fatal(FARGS, "%s", check_box(inputrec.pbcType, box));
 987     }
 988
 989     /* Test particle insertion ? */
 990     if (EI_TPI(inputrec.eI))
 991     {
 992         /* Set to the size of the molecule to be inserted (the last one) */
 993         gmx::RangePartitioning molecules = gmx_mtop_molecules(mtop);
 994         forcerec->n_tpi                  = molecules.block(molecules.numBlocks() - 1).size();
 995     }
 996     else
 997     {
 998         forcerec->n_tpi = 0;
 999     }
1000
1001     if (inputrec.coulombtype == CoulombInteractionType::RFNecUnsupported
1002         || inputrec.coulombtype == CoulombInteractionType::GRFNotused)
1003     {
1004         gmx_fatal(FARGS, "%s electrostatics is no longer supported", enumValueToString(inputrec.coulombtype));
1005     }
1006
1007     if (inputrec.bAdress)
1008     {
1009         gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1010     }
1011     if (inputrec.useTwinRange)
1012     {
1013         gmx_fatal(FARGS, "Twin-range simulations are no longer supported");
1014     }
1015     /* Copy the user determined parameters */
1016     forcerec->userint1  = inputrec.userint1;
1017     forcerec->userint2  = inputrec.userint2;
1018     forcerec->userint3  = inputrec.userint3;
1019     forcerec->userint4  = inputrec.userint4;
1020     forcerec->userreal1 = inputrec.userreal1;
1021     forcerec->userreal2 = inputrec.userreal2;
1022     forcerec->userreal3 = inputrec.userreal3;
1023     forcerec->userreal4 = inputrec.userreal4;
1024
1025     /* Shell stuff */
1026     forcerec->fc_stepsize = inputrec.fc_stepsize;
1027
1028     /* Free energy */
1029     forcerec->efep = inputrec.efep;
1030
1031     if ((getenv("GMX_DISABLE_SIMD_KERNELS") != nullptr) || (getenv("GMX_NOOPTIMIZEDKERNELS") != nullptr))
1032     {
1033         forcerec->use_simd_kernels = FALSE;
1034         if (fplog != nullptr)
1035         {
1036             fprintf(fplog,
1037                     "\nFound environment variable GMX_DISABLE_SIMD_KERNELS.\n"
1038                     "Disabling the usage of any SIMD-specific non-bonded & bonded kernel routines\n"
1039                     "(e.g. SSE2/SSE4.1/AVX).\n\n");
1040         }
1041     }
1042
1043     forcerec->haveBuckingham = (mtop.ffparams.functype[0] == F_BHAM);
1044
1045     /* Neighbour searching stuff */
1046     forcerec->pbcType = inputrec.pbcType;
1047
1048     /* Determine if we will do PBC for distances in bonded interactions */
1049     if (forcerec->pbcType == PbcType::No)
1050     {
1051         forcerec->bMolPBC = FALSE;
1052     }
1053     else
1054     {
1055         const bool useEwaldSurfaceCorrection =
1056                 (EEL_PME_EWALD(inputrec.coulombtype) && inputrec.epsilon_surface != 0);
1057         const bool haveOrientationRestraints = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
1058         if (!DOMAINDECOMP(commrec))
1059         {
1060             forcerec->bMolPBC = true;
1061
1062             if (useEwaldSurfaceCorrection || haveOrientationRestraints)
1063             {
1064                 forcerec->wholeMoleculeTransform =
1065                         std::make_unique<gmx::WholeMoleculeTransform>(mtop, inputrec.pbcType);
1066             }
1067         }
1068         else
1069         {
1070             forcerec->bMolPBC = dd_bonded_molpbc(*commrec->dd, forcerec->pbcType);
1071
1072             /* With Ewald surface correction it is useful to support e.g. running water
1073              * in parallel with update groups.
1074              * With orientation restraints there is no sensible use case supported with DD.
1075              */
1076             if ((useEwaldSurfaceCorrection && !dd_moleculesAreAlwaysWhole(*commrec->dd))
1077                 || haveOrientationRestraints)
1078             {
1079                 gmx_fatal(FARGS,
1080                           "You requested Ewald surface correction or orientation restraints, "
1081                           "but molecules are broken "
1082                           "over periodic boundary conditions by the domain decomposition. "
1083                           "Run without domain decomposition instead.");
1084             }
1085         }
1086
1087         if (useEwaldSurfaceCorrection)
1088         {
1089             GMX_RELEASE_ASSERT(!DOMAINDECOMP(commrec) || dd_moleculesAreAlwaysWhole(*commrec->dd),
1090                                "Molecules can not be broken by PBC with epsilon_surface > 0");
1091         }
1092     }
1093
1094     forcerec->rc_scaling = inputrec.refcoord_scaling;
1095     copy_rvec(inputrec.posres_com, forcerec->posres_com);
1096     copy_rvec(inputrec.posres_comB, forcerec->posres_comB);
1097     forcerec->rlist                  = cutoff_inf(inputrec.rlist);
1098     forcerec->ljpme_combination_rule = inputrec.ljpme_combination_rule;
1099
1100     /* This now calculates sum for q and c6*/
1101     bool systemHasNetCharge = set_chargesum(fplog, forcerec, mtop);
1102
1103     /* Make data structure used by kernels */
1104     forcerec->ic = std::make_unique<interaction_const_t>(
1105             init_interaction_const(fplog, inputrec, mtop, systemHasNetCharge));
1106     init_interaction_const_tables(fplog, forcerec->ic.get(), forcerec->rlist, inputrec.tabext);
1107
1108     const interaction_const_t* interactionConst = forcerec->ic.get();
1109
1110     /* TODO: Replace this Ewald table or move it into interaction_const_t */
1111     if (inputrec.coulombtype == CoulombInteractionType::Ewald)
1112     {
1113         forcerec->ewald_table = std::make_unique<gmx_ewald_tab_t>(inputrec, fplog);
1114     }
1115
1116     /* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */
1117     switch (interactionConst->eeltype)
1118     {
1119         case CoulombInteractionType::Cut:
1120             forcerec->nbkernel_elec_interaction = NbkernelElecType::Coulomb;
1121             break;
1122
1123         case CoulombInteractionType::RF:
1124         case CoulombInteractionType::RFZero:
1125             forcerec->nbkernel_elec_interaction = NbkernelElecType::ReactionField;
1126             break;
1127
1128         case CoulombInteractionType::Switch:
1129         case CoulombInteractionType::Shift:
1130         case CoulombInteractionType::User:
1131         case CoulombInteractionType::PmeSwitch:
1132         case CoulombInteractionType::PmeUser:
1133         case CoulombInteractionType::PmeUserSwitch:
1134             forcerec->nbkernel_elec_interaction = NbkernelElecType::CubicSplineTable;
1135             break;
1136
1137         case CoulombInteractionType::Pme:
1138         case CoulombInteractionType::P3mAD:
1139         case CoulombInteractionType::Ewald:
1140             forcerec->nbkernel_elec_interaction = NbkernelElecType::Ewald;
1141             break;
1142
1143         default:
1144             gmx_fatal(FARGS,
1145                       "Unsupported electrostatic interaction: %s",
1146                       enumValueToString(interactionConst->eeltype));
1147     }
1148     forcerec->nbkernel_elec_modifier = interactionConst->coulomb_modifier;
1149
1150     /* Vdw: Translate from mdp settings to kernel format */
1151     switch (interactionConst->vdwtype)
1152     {
1153         case VanDerWaalsType::Cut:
1154             if (forcerec->haveBuckingham)
1155             {
1156                 forcerec->nbkernel_vdw_interaction = NbkernelVdwType::Buckingham;
1157             }
1158             else
1159             {
1160                 forcerec->nbkernel_vdw_interaction = NbkernelVdwType::LennardJones;
1161             }
1162             break;
1163         case VanDerWaalsType::Pme:
1164             forcerec->nbkernel_vdw_interaction = NbkernelVdwType::LJEwald;
1165             break;
1166
1167         case VanDerWaalsType::Switch:
1168         case VanDerWaalsType::Shift:
1169         case VanDerWaalsType::User:
1170             forcerec->nbkernel_vdw_interaction = NbkernelVdwType::CubicSplineTable;
1171             break;
1172
1173         default:
1174             gmx_fatal(FARGS, "Unsupported vdw interaction: %s", enumValueToString(interactionConst->vdwtype));
1175     }
1176     forcerec->nbkernel_vdw_modifier = interactionConst->vdw_modifier;
1177
1178     if (!gmx_within_tol(interactionConst->reppow, 12.0, 10 * GMX_DOUBLE_EPS))
1179     {
1180         gmx_fatal(FARGS, "Only LJ repulsion power 12 is supported");
1181     }
1182     /* Older tpr files can contain Coulomb user tables with the Verlet cutoff-scheme,
1183      * while mdrun does not (and never did) support this.
1184      */
1185     if (EEL_USER(forcerec->ic->eeltype))
1186     {
1187         gmx_fatal(FARGS,
1188                   "Electrostatics type %s is currently not supported",
1189                   enumValueToString(inputrec.coulombtype));
1190     }
1191
1192     /* 1-4 interaction electrostatics */
1193     forcerec->fudgeQQ = mtop.ffparams.fudgeQQ;
1194
1195     // Multiple time stepping
1196     forcerec->useMts = inputrec.useMts;
1197
1198     if (forcerec->useMts)
1199     {
1200         GMX_ASSERT(gmx::checkMtsRequirements(inputrec).empty(),
1201                    "All MTS requirements should be met here");
1202     }
1203
1204     const bool haveDirectVirialContributionsFast =
1205             forcerec->forceProviders->hasForceProvider() || gmx_mtop_ftype_count(mtop, F_POSRES) > 0
1206             || gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0 || inputrec.nwall > 0 || inputrec.bPull
1207             || inputrec.bRot || inputrec.bIMD;
1208     const bool haveDirectVirialContributionsSlow =
1209             EEL_FULL(interactionConst->eeltype) || EVDW_PME(interactionConst->vdwtype);
1210     for (int i = 0; i < (forcerec->useMts ? 2 : 1); i++)
1211     {
1212         bool haveDirectVirialContributions =
1213                 (((!forcerec->useMts || i == 0) && haveDirectVirialContributionsFast)
1214                  || ((!forcerec->useMts || i == 1) && haveDirectVirialContributionsSlow));
1215         forcerec->forceHelperBuffers.emplace_back(haveDirectVirialContributions);
1216     }
1217
1218     if (forcerec->shift_vec.empty())
1219     {
1220         forcerec->shift_vec.resize(SHIFTS);
1221     }
1222
1223     if (forcerec->nbfp.empty())
1224     {
1225         forcerec->ntype = mtop.ffparams.atnr;
1226         forcerec->nbfp  = makeNonBondedParameterLists(mtop.ffparams, forcerec->haveBuckingham);
1227         if (EVDW_PME(interactionConst->vdwtype))
1228         {
1229             forcerec->ljpme_c6grid = makeLJPmeC6GridCorrectionParameters(mtop.ffparams, *forcerec);
1230         }
1231     }
1232
1233     /* Copy the energy group exclusions */
1234     forcerec->egp_flags = inputrec.opts.egp_flags;
1235
1236     /* Van der Waals stuff */
1237     if ((interactionConst->vdwtype != VanDerWaalsType::Cut)
1238         && (interactionConst->vdwtype != VanDerWaalsType::User) && !forcerec->haveBuckingham)
1239     {
1240         if (interactionConst->rvdw_switch >= interactionConst->rvdw)
1241         {
1242             gmx_fatal(FARGS,
1243                       "rvdw_switch (%f) must be < rvdw (%f)",
1244                       interactionConst->rvdw_switch,
1245                       interactionConst->rvdw);
1246         }
1247         if (fplog)
1248         {
1249             fprintf(fplog,
1250                     "Using %s Lennard-Jones, switch between %g and %g nm\n",
1251                     (interactionConst->eeltype == CoulombInteractionType::Switch) ? "switched" : "shifted",
1252                     interactionConst->rvdw_switch,
1253                     interactionConst->rvdw);
1254         }
1255     }
1256
1257     if (forcerec->haveBuckingham && EVDW_PME(interactionConst->vdwtype))
1258     {
1259         gmx_fatal(FARGS, "LJ PME not supported with Buckingham");
1260     }
1261
1262     if (forcerec->haveBuckingham
1263         && (interactionConst->vdwtype == VanDerWaalsType::Shift
1264             || interactionConst->vdwtype == VanDerWaalsType::Switch))
1265     {
1266         gmx_fatal(FARGS, "Switch/shift interaction not supported with Buckingham");
1267     }
1268
1269     if (forcerec->haveBuckingham)
1270     {
1271         gmx_fatal(FARGS, "The Verlet cutoff-scheme does not (yet) support Buckingham");
1272     }
1273
1274     if (inputrec.implicit_solvent)
1275     {
1276         gmx_fatal(FARGS, "Implict solvation is no longer supported.");
1277     }
1278
1279
1280     /* This code automatically gives table length tabext without cut-off's,
1281      * in that case grompp should already have checked that we do not need
1282      * normal tables and we only generate tables for 1-4 interactions.
1283      */
1284     real rtab = inputrec.rlist + inputrec.tabext;
1285
1286     /* We want to use unmodified tables for 1-4 coulombic
1287      * interactions, so we must in general have an extra set of
1288      * tables. */
1289     if (gmx_mtop_ftype_count(mtop, F_LJ14) > 0 || gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0
1290         || gmx_mtop_ftype_count(mtop, F_LJC_PAIRS_NB) > 0)
1291     {
1292         forcerec->pairsTable = make_tables(fplog, interactionConst, tabpfn, rtab, GMX_MAKETABLES_14ONLY);
1293     }
1294
1295     /* Wall stuff */
1296     forcerec->nwall = inputrec.nwall;
1297     if (inputrec.nwall && inputrec.wall_type == WallType::Table)
1298     {
1299         make_wall_tables(fplog, inputrec, tabfn, &mtop.groups, forcerec);
1300     }
1301
1302     forcerec->fcdata = std::make_unique<t_fcdata>();
1303
1304     if (!tabbfnm.empty())
1305     {
1306         t_fcdata& fcdata = *forcerec->fcdata;
1307         // Need to catch std::bad_alloc
1308         // TODO Don't need to catch this here, when merging with master branch
1309         try
1310         {
1311             // TODO move these tables into a separate struct and store reference in ListedForces
1312             fcdata.bondtab = make_bonded_tables(fplog, F_TABBONDS, F_TABBONDSNC, mtop, tabbfnm, "b");
1313             fcdata.angletab = make_bonded_tables(fplog, F_TABANGLES, -1, mtop, tabbfnm, "a");
1314             fcdata.dihtab   = make_bonded_tables(fplog, F_TABDIHS, -1, mtop, tabbfnm, "d");
1315         }
1316         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1317     }
1318     else
1319     {
1320         if (debug)
1321         {
1322             fprintf(debug,
1323                     "No fcdata or table file name passed, can not read table, can not do bonded "
1324                     "interactions\n");
1325         }
1326     }
1327
1328     /* Initialize the thread working data for bonded interactions */
1329     if (forcerec->useMts)
1330     {
1331         // Add one ListedForces object for each MTS level
1332         bool isFirstLevel = true;
1333         for (const auto& mtsLevel : inputrec.mtsLevels)
1334         {
1335             ListedForces::InteractionSelection interactionSelection;
1336             const auto&                        forceGroups = mtsLevel.forceGroups;
1337             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Pair)])
1338             {
1339                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Pairs));
1340             }
1341             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Dihedral)])
1342             {
1343                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Dihedrals));
1344             }
1345             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Angle)])
1346             {
1347                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Angles));
1348             }
1349             if (isFirstLevel)
1350             {
1351                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Rest));
1352                 isFirstLevel = false;
1353             }
1354             forcerec->listedForces.emplace_back(
1355                     mtop.ffparams,
1356                     mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1357                     gmx_omp_nthreads_get(emntBonded),
1358                     interactionSelection,
1359                     fplog);
1360         }
1361     }
1362     else
1363     {
1364         // Add one ListedForces object with all listed interactions
1365         forcerec->listedForces.emplace_back(
1366                 mtop.ffparams,
1367                 mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1368                 gmx_omp_nthreads_get(emntBonded),
1369                 ListedForces::interactionSelectionAll(),
1370                 fplog);
1371     }
1372
1373     // QM/MM initialization if requested
1374     if (inputrec.bQMMM)
1375     {
1376         gmx_incons("QM/MM was requested, but is no longer available in GROMACS");
1377     }
1378
1379     /* Set all the static charge group info */
1380     forcerec->cginfo_mb = init_cginfo_mb(mtop, forcerec);
1381     if (!DOMAINDECOMP(commrec))
1382     {
1383         forcerec->cginfo = cginfo_expand(mtop.molblock.size(), forcerec->cginfo_mb);
1384     }
1385
1386     if (!DOMAINDECOMP(commrec))
1387     {
1388         forcerec_set_ranges(forcerec, mtop.natoms, mtop.natoms, mtop.natoms);
1389     }
1390
1391     forcerec->print_force = print_force;
1392
1393     forcerec->nthread_ewc = gmx_omp_nthreads_get(emntBonded);
1394     forcerec->ewc_t.resize(forcerec->nthread_ewc);
1395
1396     if (inputrec.eDispCorr != DispersionCorrectionType::No)
1397     {
1398         forcerec->dispersionCorrection = std::make_unique<DispersionCorrection>(
1399                 mtop, inputrec, forcerec->haveBuckingham, forcerec->ntype, forcerec->nbfp, *forcerec->ic, tabfn);
1400         forcerec->dispersionCorrection->print(mdlog);
1401     }
1402
1403     if (fplog != nullptr)
1404     {
1405         /* Here we switch from using mdlog, which prints the newline before
1406          * the paragraph, to our old fprintf logging, which prints the newline
1407          * after the paragraph, so we should add a newline here.
1408          */
1409         fprintf(fplog, "\n");
1410     }
1411 }
1412
1413 t_forcerec::t_forcerec() = default;
1414
1415 t_forcerec::~t_forcerec() = default;