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 <cmath>
  44 #include <cstdlib>
  45 #include <cstring>
  46
  47 #include <algorithm>
  48 #include <memory>
  49
  50 #include "gromacs/commandline/filenm.h"
  51 #include "gromacs/domdec/domdec.h"
  52 #include "gromacs/domdec/domdec_struct.h"
  53 #include "gromacs/ewald/ewald.h"
  54 #include "gromacs/ewald/pme_pp_comm_gpu.h"
  55 #include "gromacs/fileio/filetypes.h"
  56 #include "gromacs/gmxlib/network.h"
  57 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  58 #include "gromacs/gpu_utils/gpu_utils.h"
  59 #include "gromacs/hardware/hw_info.h"
  60 #include "gromacs/listed_forces/gpubonded.h"
  61 #include "gromacs/listed_forces/listed_forces.h"
  62 #include "gromacs/listed_forces/pairs.h"
  63 #include "gromacs/math/functions.h"
  64 #include "gromacs/math/utilities.h"
  65 #include "gromacs/math/vec.h"
  66 #include "gromacs/mdlib/dispersioncorrection.h"
  67 #include "gromacs/mdlib/force.h"
  68 #include "gromacs/mdlib/forcerec_threading.h"
  69 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  70 #include "gromacs/mdlib/md_support.h"
  71 #include "gromacs/mdlib/wall.h"
  72 #include "gromacs/mdlib/wholemoleculetransform.h"
  73 #include "gromacs/mdtypes/commrec.h"
  74 #include "gromacs/mdtypes/fcdata.h"
  75 #include "gromacs/mdtypes/forcerec.h"
  76 #include "gromacs/mdtypes/group.h"
  77 #include "gromacs/mdtypes/iforceprovider.h"
  78 #include "gromacs/mdtypes/inputrec.h"
  79 #include "gromacs/mdtypes/interaction_const.h"
  80 #include "gromacs/mdtypes/md_enums.h"
  81 #include "gromacs/mdtypes/multipletimestepping.h"
  82 #include "gromacs/mdtypes/nblist.h"
  83 #include "gromacs/nbnxm/nbnxm.h"
  84 #include "gromacs/pbcutil/ishift.h"
  85 #include "gromacs/pbcutil/pbc.h"
  86 #include "gromacs/tables/forcetable.h"
  87 #include "gromacs/topology/mtop_util.h"
  88 #include "gromacs/topology/idef.h"
  89 #include "gromacs/trajectory/trajectoryframe.h"
  90 #include "gromacs/utility/cstringutil.h"
  91 #include "gromacs/utility/exceptions.h"
  92 #include "gromacs/utility/fatalerror.h"
  93 #include "gromacs/utility/gmxassert.h"
  94 #include "gromacs/utility/logger.h"
  95 #include "gromacs/utility/physicalnodecommunicator.h"
  96 #include "gromacs/utility/pleasecite.h"
  97 #include "gromacs/utility/smalloc.h"
  98 #include "gromacs/utility/strconvert.h"
  99
 100 #include "gpuforcereduction.h"
 101
 102 ForceHelperBuffers::ForceHelperBuffers(bool haveDirectVirialContributions) :
 103     haveDirectVirialContributions_(haveDirectVirialContributions)
 104 {
 105     shiftForces_.resize(gmx::c_numShiftVectors);
 106 }
 107
 108 void ForceHelperBuffers::resize(int numAtoms)
 109 {
 110     if (haveDirectVirialContributions_)
 111     {
 112         forceBufferForDirectVirialContributions_.resize(numAtoms);
 113     }
 114 }
 115
 116 std::vector<real> makeNonBondedParameterLists(const int                      numAtomTypes,
 117                                               gmx::ArrayRef<const t_iparams> iparams,
 118                                               bool                           useBuckinghamPotential)
 119 {
 120     std::vector<real> nbfp;
 121
 122     if (useBuckinghamPotential)
 123     {
 124         nbfp.resize(3 * numAtomTypes * numAtomTypes);
 125         int k = 0;
 126         for (int i = 0; (i < numAtomTypes); i++)
 127         {
 128             for (int j = 0; (j < numAtomTypes); j++, k++)
 129             {
 130                 BHAMA(nbfp, numAtomTypes, i, j) = iparams[k].bham.a;
 131                 BHAMB(nbfp, numAtomTypes, i, j) = iparams[k].bham.b;
 132                 /* nbfp now includes the 6.0 derivative prefactor */
 133                 BHAMC(nbfp, numAtomTypes, i, j) = iparams[k].bham.c * 6.0;
 134             }
 135         }
 136     }
 137     else
 138     {
 139         nbfp.resize(2 * numAtomTypes * numAtomTypes);
 140         int k = 0;
 141         for (int i = 0; (i < numAtomTypes); i++)
 142         {
 143             for (int j = 0; (j < numAtomTypes); j++, k++)
 144             {
 145                 /* nbfp now includes the 6.0/12.0 derivative prefactors */
 146                 C6(nbfp, numAtomTypes, i, j)  = iparams[k].lj.c6 * 6.0;
 147                 C12(nbfp, numAtomTypes, i, j) = iparams[k].lj.c12 * 12.0;
 148             }
 149         }
 150     }
 151
 152     return nbfp;
 153 }
 154
 155 std::vector<real> makeLJPmeC6GridCorrectionParameters(const int                      numAtomTypes,
 156                                                       gmx::ArrayRef<const t_iparams> iparams,
 157                                                       LongRangeVdW ljpme_combination_rule)
 158 {
 159     /* For LJ-PME simulations, we correct the energies with the reciprocal space
 160      * inside of the cut-off. To do this the non-bonded kernels needs to have
 161      * access to the C6-values used on the reciprocal grid in pme.c
 162      */
 163
 164     std::vector<real> grid(2 * numAtomTypes * numAtomTypes, 0.0);
 165     int               k = 0;
 166     for (int i = 0; (i < numAtomTypes); i++)
 167     {
 168         for (int j = 0; (j < numAtomTypes); j++, k++)
 169         {
 170             real c6i  = iparams[i * (numAtomTypes + 1)].lj.c6;
 171             real c12i = iparams[i * (numAtomTypes + 1)].lj.c12;
 172             real c6j  = iparams[j * (numAtomTypes + 1)].lj.c6;
 173             real c12j = iparams[j * (numAtomTypes + 1)].lj.c12;
 174             real c6   = std::sqrt(c6i * c6j);
 175             if (ljpme_combination_rule == LongRangeVdW::LB && !gmx_numzero(c6) && !gmx_numzero(c12i)
 176                 && !gmx_numzero(c12j))
 177             {
 178                 real sigmai = gmx::sixthroot(c12i / c6i);
 179                 real sigmaj = gmx::sixthroot(c12j / c6j);
 180                 real epsi   = c6i * c6i / c12i;
 181                 real epsj   = c6j * c6j / c12j;
 182                 c6          = std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
 183             }
 184             /* Store the elements at the same relative positions as C6 in nbfp in order
 185              * to simplify access in the kernels
 186              */
 187             grid[2 * (numAtomTypes * i + j)] = c6 * 6.0;
 188         }
 189     }
 190     return grid;
 191 }
 192
 193 enum
 194 {
 195     acNONE = 0,
 196     acCONSTRAINT,
 197     acSETTLE
 198 };
 199
 200 static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t& mtop, const t_forcerec* fr)
 201 {
 202     gmx_bool* type_VDW;
 203     int*      a_con;
 204
 205     snew(type_VDW, fr->ntype);
 206     for (int ai = 0; ai < fr->ntype; ai++)
 207     {
 208         type_VDW[ai] = FALSE;
 209         for (int j = 0; j < fr->ntype; j++)
 210         {
 211             type_VDW[ai] = type_VDW[ai] || fr->haveBuckingham || C6(fr->nbfp, fr->ntype, ai, j) != 0
 212                            || C12(fr->nbfp, fr->ntype, ai, j) != 0;
 213         }
 214     }
 215
 216     std::vector<cginfo_mb_t> cginfoPerMolblock;
 217     int                      a_offset = 0;
 218     for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
 219     {
 220         const gmx_molblock_t& molb = mtop.molblock[mb];
 221         const gmx_moltype_t&  molt = mtop.moltype[molb.type];
 222         const auto&           excl = molt.excls;
 223
 224         /* Check if the cginfo is identical for all molecules in this block.
 225          * If so, we only need an array of the size of one molecule.
 226          * Otherwise we make an array of #mol times #cgs per molecule.
 227          */
 228         gmx_bool bId = TRUE;
 229         for (int m = 0; m < molb.nmol; m++)
 230         {
 231             const int am = m * molt.atoms.nr;
 232             for (int a = 0; a < molt.atoms.nr; a++)
 233             {
 234                 if (getGroupType(mtop.groups, SimulationAtomGroupType::QuantumMechanics, a_offset + am + a)
 235                     != getGroupType(mtop.groups, SimulationAtomGroupType::QuantumMechanics, a_offset + a))
 236                 {
 237                     bId = FALSE;
 238                 }
 239                 if (!mtop.groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].empty())
 240                 {
 241                     if (mtop.groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics][a_offset + am + a]
 242                         != mtop.groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics][a_offset + a])
 243                     {
 244                         bId = FALSE;
 245                     }
 246                 }
 247             }
 248         }
 249
 250         cginfo_mb_t cginfo_mb;
 251         cginfo_mb.cg_start = a_offset;
 252         cginfo_mb.cg_end   = a_offset + molb.nmol * molt.atoms.nr;
 253         cginfo_mb.cg_mod   = (bId ? 1 : molb.nmol) * molt.atoms.nr;
 254         cginfo_mb.cginfo.resize(cginfo_mb.cg_mod);
 255         gmx::ArrayRef<int> cginfo = cginfo_mb.cginfo;
 256
 257         /* Set constraints flags for constrained atoms */
 258         snew(a_con, molt.atoms.nr);
 259         for (int ftype = 0; ftype < F_NRE; ftype++)
 260         {
 261             if (interaction_function[ftype].flags & IF_CONSTRAINT)
 262             {
 263                 const int nral = NRAL(ftype);
 264                 for (int ia = 0; ia < molt.ilist[ftype].size(); ia += 1 + nral)
 265                 {
 266                     int a;
 267
 268                     for (a = 0; a < nral; a++)
 269                     {
 270                         a_con[molt.ilist[ftype].iatoms[ia + 1 + a]] =
 271                                 (ftype == F_SETTLE ? acSETTLE : acCONSTRAINT);
 272                     }
 273                 }
 274             }
 275         }
 276
 277         for (int m = 0; m < (bId ? 1 : molb.nmol); m++)
 278         {
 279             const int molculeOffsetInBlock = m * molt.atoms.nr;
 280             for (int a = 0; a < molt.atoms.nr; a++)
 281             {
 282                 const t_atom& atom     = molt.atoms.atom[a];
 283                 int&          atomInfo = cginfo[molculeOffsetInBlock + a];
 284
 285                 /* Store the energy group in cginfo */
 286                 int gid = getGroupType(mtop.groups,
 287                                        SimulationAtomGroupType::EnergyOutput,
 288                                        a_offset + molculeOffsetInBlock + a);
 289                 SET_CGINFO_GID(atomInfo, gid);
 290
 291                 bool bHaveVDW = (type_VDW[atom.type] || type_VDW[atom.typeB]);
 292                 bool bHaveQ   = (atom.q != 0 || atom.qB != 0);
 293
 294                 bool haveExclusions = false;
 295                 /* Loop over all the exclusions of atom ai */
 296                 for (const int j : excl[a])
 297                 {
 298                     if (j != a)
 299                     {
 300                         haveExclusions = true;
 301                         break;
 302                     }
 303                 }
 304
 305                 switch (a_con[a])
 306                 {
 307                     case acCONSTRAINT: SET_CGINFO_CONSTR(atomInfo); break;
 308                     case acSETTLE: SET_CGINFO_SETTLE(atomInfo); break;
 309                     default: break;
 310                 }
 311                 if (haveExclusions)
 312                 {
 313                     SET_CGINFO_EXCL_INTER(atomInfo);
 314                 }
 315                 if (bHaveVDW)
 316                 {
 317                     SET_CGINFO_HAS_VDW(atomInfo);
 318                 }
 319                 if (bHaveQ)
 320                 {
 321                     SET_CGINFO_HAS_Q(atomInfo);
 322                 }
 323                 if (fr->efep != FreeEnergyPerturbationType::No && PERTURBED(atom))
 324                 {
 325                     SET_CGINFO_FEP(atomInfo);
 326                 }
 327             }
 328         }
 329
 330         sfree(a_con);
 331
 332         cginfoPerMolblock.push_back(cginfo_mb);
 333
 334         a_offset += molb.nmol * molt.atoms.nr;
 335     }
 336     sfree(type_VDW);
 337
 338     return cginfoPerMolblock;
 339 }
 340
 341 static std::vector<int> cginfo_expand(const int nmb, gmx::ArrayRef<const cginfo_mb_t> cgi_mb)
 342 {
 343     const int ncg = cgi_mb[nmb - 1].cg_end;
 344
 345     std::vector<int> cginfo(ncg);
 346
 347     int mb = 0;
 348     for (int cg = 0; cg < ncg; cg++)
 349     {
 350         while (cg >= cgi_mb[mb].cg_end)
 351         {
 352             mb++;
 353         }
 354         cginfo[cg] = cgi_mb[mb].cginfo[(cg - cgi_mb[mb].cg_start) % cgi_mb[mb].cg_mod];
 355     }
 356
 357     return cginfo;
 358 }
 359
 360 /* Sets the sum of charges (squared) and C6 in the system in fr.
 361  * Returns whether the system has a net charge.
 362  */
 363 static bool set_chargesum(FILE* log, t_forcerec* fr, const gmx_mtop_t& mtop)
 364 {
 365     /*This now calculates sum for q and c6*/
 366     double qsum, q2sum, q, c6sum, c6;
 367
 368     qsum  = 0;
 369     q2sum = 0;
 370     c6sum = 0;
 371     for (const gmx_molblock_t& molb : mtop.molblock)
 372     {
 373         int            nmol  = molb.nmol;
 374         const t_atoms* atoms = &mtop.moltype[molb.type].atoms;
 375         for (int i = 0; i < atoms->nr; i++)
 376         {
 377             q = atoms->atom[i].q;
 378             qsum += nmol * q;
 379             q2sum += nmol * q * q;
 380             c6 = mtop.ffparams.iparams[atoms->atom[i].type * (mtop.ffparams.atnr + 1)].lj.c6;
 381             c6sum += nmol * c6;
 382         }
 383     }
 384     fr->qsum[0]  = qsum;
 385     fr->q2sum[0] = q2sum;
 386     fr->c6sum[0] = c6sum;
 387
 388     if (fr->efep != FreeEnergyPerturbationType::No)
 389     {
 390         qsum  = 0;
 391         q2sum = 0;
 392         c6sum = 0;
 393         for (const gmx_molblock_t& molb : mtop.molblock)
 394         {
 395             int            nmol  = molb.nmol;
 396             const t_atoms* atoms = &mtop.moltype[molb.type].atoms;
 397             for (int i = 0; i < atoms->nr; i++)
 398             {
 399                 q = atoms->atom[i].qB;
 400                 qsum += nmol * q;
 401                 q2sum += nmol * q * q;
 402                 c6 = mtop.ffparams.iparams[atoms->atom[i].typeB * (mtop.ffparams.atnr + 1)].lj.c6;
 403                 c6sum += nmol * c6;
 404             }
 405             fr->qsum[1]  = qsum;
 406             fr->q2sum[1] = q2sum;
 407             fr->c6sum[1] = c6sum;
 408         }
 409     }
 410     else
 411     {
 412         fr->qsum[1]  = fr->qsum[0];
 413         fr->q2sum[1] = fr->q2sum[0];
 414         fr->c6sum[1] = fr->c6sum[0];
 415     }
 416     if (log)
 417     {
 418         if (fr->efep == FreeEnergyPerturbationType::No)
 419         {
 420             fprintf(log, "System total charge: %.3f\n", fr->qsum[0]);
 421         }
 422         else
 423         {
 424             fprintf(log, "System total charge, top. A: %.3f top. B: %.3f\n", fr->qsum[0], fr->qsum[1]);
 425         }
 426     }
 427
 428     /* A cut-off of 1e-4 is used to catch rounding errors due to ascii input */
 429     return (std::abs(fr->qsum[0]) > 1e-4 || std::abs(fr->qsum[1]) > 1e-4);
 430 }
 431
 432 /*!\brief If there's bonded interactions of type \c ftype1 or \c
 433  * ftype2 present in the topology, build an array of the number of
 434  * interactions present for each bonded interaction index found in the
 435  * topology.
 436  *
 437  * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
 438  * valid type with that parameter.
 439  *
 440  * \c count will be reallocated as necessary to fit the largest bonded
 441  * interaction index found, and its current size will be returned in
 442  * \c ncount. It will contain zero for every bonded interaction index
 443  * for which no interactions are present in the topology.
 444  */
 445 static void count_tables(int ftype1, int ftype2, const gmx_mtop_t& mtop, int* ncount, int** count)
 446 {
 447     int ftype, i, j, tabnr;
 448
 449     // Loop over all moleculetypes
 450     for (const gmx_moltype_t& molt : mtop.moltype)
 451     {
 452         // Loop over all interaction types
 453         for (ftype = 0; ftype < F_NRE; ftype++)
 454         {
 455             // If the current interaction type is one of the types whose tables we're trying to count...
 456             if (ftype == ftype1 || ftype == ftype2)
 457             {
 458                 const InteractionList& il     = molt.ilist[ftype];
 459                 const int              stride = 1 + NRAL(ftype);
 460                 // ... and there are actually some interactions for this type
 461                 for (i = 0; i < il.size(); i += stride)
 462                 {
 463                     // Find out which table index the user wanted
 464                     tabnr = mtop.ffparams.iparams[il.iatoms[i]].tab.table;
 465                     if (tabnr < 0)
 466                     {
 467                         gmx_fatal(FARGS, "A bonded table number is smaller than 0: %d\n", tabnr);
 468                     }
 469                     // Make room for this index in the data structure
 470                     if (tabnr >= *ncount)
 471                     {
 472                         srenew(*count, tabnr + 1);
 473                         for (j = *ncount; j < tabnr + 1; j++)
 474                         {
 475                             (*count)[j] = 0;
 476                         }
 477                         *ncount = tabnr + 1;
 478                     }
 479                     // Record that this table index is used and must have a valid file
 480                     (*count)[tabnr]++;
 481                 }
 482             }
 483         }
 484     }
 485 }
 486
 487 /*!\brief If there's bonded interactions of flavour \c tabext and type
 488  * \c ftype1 or \c ftype2 present in the topology, seek them in the
 489  * list of filenames passed to mdrun, and make bonded tables from
 490  * those files.
 491  *
 492  * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
 493  * valid type with that parameter.
 494  *
 495  * A fatal error occurs if no matching filename is found.
 496  */
 497 static std::vector<bondedtable_t> make_bonded_tables(FILE*                            fplog,
 498                                                      int                              ftype1,
 499                                                      int                              ftype2,
 500                                                      const gmx_mtop_t&                mtop,
 501                                                      gmx::ArrayRef<const std::string> tabbfnm,
 502                                                      const char*                      tabext)
 503 {
 504     std::vector<bondedtable_t> tab;
 505
 506     int  ncount = 0;
 507     int* count  = nullptr;
 508     count_tables(ftype1, ftype2, mtop, &ncount, &count);
 509
 510     // Are there any relevant tabulated bond interactions?
 511     if (ncount > 0)
 512     {
 513         tab.resize(ncount);
 514         for (int i = 0; i < ncount; i++)
 515         {
 516             // Do any interactions exist that requires this table?
 517             if (count[i] > 0)
 518             {
 519                 // This pattern enforces the current requirement that
 520                 // table filenames end in a characteristic sequence
 521                 // before the file type extension, and avoids table 13
 522                 // being recognized and used for table 1.
 523                 std::string patternToFind = gmx::formatString("_%s%d.%s", tabext, i, ftp2ext(efXVG));
 524                 bool        madeTable     = false;
 525                 for (gmx::index j = 0; j < tabbfnm.ssize() && !madeTable; ++j)
 526                 {
 527                     if (gmx::endsWith(tabbfnm[j], patternToFind))
 528                     {
 529                         // Finally read the table from the file found
 530                         tab[i]    = make_bonded_table(fplog, tabbfnm[j].c_str(), NRAL(ftype1) - 2);
 531                         madeTable = true;
 532                     }
 533                 }
 534                 if (!madeTable)
 535                 {
 536                     bool isPlural = (ftype2 != -1);
 537                     gmx_fatal(FARGS,
 538                               "Tabulated interaction of type '%s%s%s' with index %d cannot be used "
 539                               "because no table file whose name matched '%s' was passed via the "
 540                               "gmx mdrun -tableb command-line option.",
 541                               interaction_function[ftype1].longname,
 542                               isPlural ? "' or '" : "",
 543                               isPlural ? interaction_function[ftype2].longname : "",
 544                               i,
 545                               patternToFind.c_str());
 546                 }
 547             }
 548         }
 549         sfree(count);
 550     }
 551
 552     return tab;
 553 }
 554
 555 void forcerec_set_ranges(t_forcerec* fr, int natoms_force, int natoms_force_constr, int natoms_f_novirsum)
 556 {
 557     fr->natoms_force        = natoms_force;
 558     fr->natoms_force_constr = natoms_force_constr;
 559
 560     for (auto& forceHelperBuffers : fr->forceHelperBuffers)
 561     {
 562         forceHelperBuffers.resize(natoms_f_novirsum);
 563     }
 564 }
 565
 566 /* Generate Coulomb and/or Van der Waals Ewald long-range correction tables
 567  *
 568  * Tables are generated for one or both, depending on if the pointers
 569  * are non-null. The spacing for both table sets is the same and obeys
 570  * both accuracy requirements, when relevant.
 571  */
 572 static void init_ewald_f_table(const interaction_const_t& ic,
 573                                const real                 rlist,
 574                                const real                 tabext,
 575                                EwaldCorrectionTables*     coulombTables,
 576                                EwaldCorrectionTables*     vdwTables)
 577 {
 578     const bool useCoulombTable = (EEL_PME_EWALD(ic.eeltype) && coulombTables != nullptr);
 579     const bool useVdwTable     = (EVDW_PME(ic.vdwtype) && vdwTables != nullptr);
 580
 581     /* Get the Ewald table spacing based on Coulomb and/or LJ
 582      * Ewald coefficients and rtol.
 583      */
 584     const real tableScale = ewald_spline3_table_scale(ic, useCoulombTable, useVdwTable);
 585
 586     const bool havePerturbedNonbondeds = (ic.softCoreParameters != nullptr);
 587
 588     real tableLen = ic.rcoulomb;
 589     if ((useCoulombTable || useVdwTable) && havePerturbedNonbondeds && rlist + tabext > 0.0)
 590     {
 591         /* TODO: Ideally this should also check if couple-intramol == no, but that isn't
 592          * stored in ir. Grompp puts that info into an opts structure that doesn't make it into the tpr.
 593          * The alternative is to look through all the exclusions and check if they come from
 594          * couple-intramol == no. Meanwhile, always having larger tables should only affect
 595          * memory consumption, not speed (barring cache issues).
 596          */
 597         tableLen = rlist + tabext;
 598     }
 599     const int tableSize = static_cast<int>(tableLen * tableScale) + 2;
 600
 601     if (useCoulombTable)
 602     {
 603         *coulombTables =
 604                 generateEwaldCorrectionTables(tableSize, tableScale, ic.ewaldcoeff_q, v_q_ewald_lr);
 605     }
 606
 607     if (useVdwTable)
 608     {
 609         *vdwTables = generateEwaldCorrectionTables(tableSize, tableScale, ic.ewaldcoeff_lj, v_lj_ewald_lr);
 610     }
 611 }
 612
 613 void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real rlist, const real tableExtensionLength)
 614 {
 615     if (EEL_PME_EWALD(ic->eeltype) || EVDW_PME(ic->vdwtype))
 616     {
 617         init_ewald_f_table(
 618                 *ic, rlist, tableExtensionLength, ic->coulombEwaldTables.get(), ic->vdwEwaldTables.get());
 619         if (fp != nullptr)
 620         {
 621             fprintf(fp,
 622                     "Initialized non-bonded Ewald tables, spacing: %.2e size: %zu\n\n",
 623                     1 / ic->coulombEwaldTables->scale,
 624                     ic->coulombEwaldTables->tableF.size());
 625         }
 626     }
 627 }
 628
 629 real cutoff_inf(real cutoff)
 630 {
 631     if (cutoff == 0)
 632     {
 633         cutoff = GMX_CUTOFF_INF;
 634     }
 635
 636     return cutoff;
 637 }
 638
 639 void init_forcerec(FILE*                            fplog,
 640                    const gmx::MDLogger&             mdlog,
 641                    t_forcerec*                      forcerec,
 642                    const t_inputrec&                inputrec,
 643                    const gmx_mtop_t&                mtop,
 644                    const t_commrec*                 commrec,
 645                    matrix                           box,
 646                    const char*                      tabfn,
 647                    const char*                      tabpfn,
 648                    gmx::ArrayRef<const std::string> tabbfnm,
 649                    real                             print_force)
 650 {
 651     /* The CMake default turns SIMD kernels on, but it might be turned off further down... */
 652     forcerec->use_simd_kernels = GMX_USE_SIMD_KERNELS;
 653
 654     if (check_box(inputrec.pbcType, box))
 655     {
 656         gmx_fatal(FARGS, "%s", check_box(inputrec.pbcType, box));
 657     }
 658
 659     /* Test particle insertion ? */
 660     if (EI_TPI(inputrec.eI))
 661     {
 662         /* Set to the size of the molecule to be inserted (the last one) */
 663         gmx::RangePartitioning molecules = gmx_mtop_molecules(mtop);
 664         forcerec->n_tpi                  = molecules.block(molecules.numBlocks() - 1).size();
 665     }
 666     else
 667     {
 668         forcerec->n_tpi = 0;
 669     }
 670
 671     if (inputrec.coulombtype == CoulombInteractionType::RFNecUnsupported
 672         || inputrec.coulombtype == CoulombInteractionType::GRFNotused)
 673     {
 674         gmx_fatal(FARGS, "%s electrostatics is no longer supported", enumValueToString(inputrec.coulombtype));
 675     }
 676
 677     if (inputrec.bAdress)
 678     {
 679         gmx_fatal(FARGS, "AdResS simulations are no longer supported");
 680     }
 681     if (inputrec.useTwinRange)
 682     {
 683         gmx_fatal(FARGS, "Twin-range simulations are no longer supported");
 684     }
 685     /* Copy the user determined parameters */
 686     forcerec->userint1  = inputrec.userint1;
 687     forcerec->userint2  = inputrec.userint2;
 688     forcerec->userint3  = inputrec.userint3;
 689     forcerec->userint4  = inputrec.userint4;
 690     forcerec->userreal1 = inputrec.userreal1;
 691     forcerec->userreal2 = inputrec.userreal2;
 692     forcerec->userreal3 = inputrec.userreal3;
 693     forcerec->userreal4 = inputrec.userreal4;
 694
 695     /* Shell stuff */
 696     forcerec->fc_stepsize = inputrec.fc_stepsize;
 697
 698     /* Free energy */
 699     forcerec->efep = inputrec.efep;
 700
 701     if ((getenv("GMX_DISABLE_SIMD_KERNELS") != nullptr) || (getenv("GMX_NOOPTIMIZEDKERNELS") != nullptr))
 702     {
 703         forcerec->use_simd_kernels = FALSE;
 704         if (fplog != nullptr)
 705         {
 706             fprintf(fplog,
 707                     "\nFound environment variable GMX_DISABLE_SIMD_KERNELS.\n"
 708                     "Disabling the usage of any SIMD-specific non-bonded & bonded kernel routines\n"
 709                     "(e.g. SSE2/SSE4.1/AVX).\n\n");
 710         }
 711     }
 712
 713     forcerec->haveBuckingham = (mtop.ffparams.functype[0] == F_BHAM);
 714
 715     /* Neighbour searching stuff */
 716     forcerec->pbcType = inputrec.pbcType;
 717
 718     /* Determine if we will do PBC for distances in bonded interactions */
 719     if (forcerec->pbcType == PbcType::No)
 720     {
 721         forcerec->bMolPBC = FALSE;
 722     }
 723     else
 724     {
 725         const bool useEwaldSurfaceCorrection =
 726                 (EEL_PME_EWALD(inputrec.coulombtype) && inputrec.epsilon_surface != 0);
 727         const bool haveOrientationRestraints = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
 728         if (!DOMAINDECOMP(commrec))
 729         {
 730             forcerec->bMolPBC = true;
 731
 732             if (useEwaldSurfaceCorrection || haveOrientationRestraints)
 733             {
 734                 forcerec->wholeMoleculeTransform =
 735                         std::make_unique<gmx::WholeMoleculeTransform>(mtop, inputrec.pbcType);
 736             }
 737         }
 738         else
 739         {
 740             forcerec->bMolPBC = dd_bonded_molpbc(*commrec->dd, forcerec->pbcType);
 741
 742             /* With Ewald surface correction it is useful to support e.g. running water
 743              * in parallel with update groups.
 744              * With orientation restraints there is no sensible use case supported with DD.
 745              */
 746             if ((useEwaldSurfaceCorrection && !dd_moleculesAreAlwaysWhole(*commrec->dd))
 747                 || haveOrientationRestraints)
 748             {
 749                 gmx_fatal(FARGS,
 750                           "You requested Ewald surface correction or orientation restraints, "
 751                           "but molecules are broken "
 752                           "over periodic boundary conditions by the domain decomposition. "
 753                           "Run without domain decomposition instead.");
 754             }
 755         }
 756
 757         if (useEwaldSurfaceCorrection)
 758         {
 759             GMX_RELEASE_ASSERT(!DOMAINDECOMP(commrec) || dd_moleculesAreAlwaysWhole(*commrec->dd),
 760                                "Molecules can not be broken by PBC with epsilon_surface > 0");
 761         }
 762     }
 763
 764     forcerec->rc_scaling = inputrec.refcoord_scaling;
 765     copy_rvec(inputrec.posres_com, forcerec->posres_com);
 766     copy_rvec(inputrec.posres_comB, forcerec->posres_comB);
 767     forcerec->rlist                  = cutoff_inf(inputrec.rlist);
 768     forcerec->ljpme_combination_rule = inputrec.ljpme_combination_rule;
 769
 770     /* This now calculates sum for q and c6*/
 771     bool systemHasNetCharge = set_chargesum(fplog, forcerec, mtop);
 772
 773     /* Make data structure used by kernels */
 774     forcerec->ic = std::make_unique<interaction_const_t>(
 775             init_interaction_const(fplog, inputrec, mtop, systemHasNetCharge));
 776     init_interaction_const_tables(fplog, forcerec->ic.get(), forcerec->rlist, inputrec.tabext);
 777
 778     const interaction_const_t* interactionConst = forcerec->ic.get();
 779
 780     /* TODO: Replace this Ewald table or move it into interaction_const_t */
 781     if (inputrec.coulombtype == CoulombInteractionType::Ewald)
 782     {
 783         forcerec->ewald_table = std::make_unique<gmx_ewald_tab_t>(inputrec, fplog);
 784     }
 785
 786     /* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */
 787     switch (interactionConst->eeltype)
 788     {
 789         case CoulombInteractionType::Cut:
 790             forcerec->nbkernel_elec_interaction = NbkernelElecType::Coulomb;
 791             break;
 792
 793         case CoulombInteractionType::RF:
 794         case CoulombInteractionType::RFZero:
 795             forcerec->nbkernel_elec_interaction = NbkernelElecType::ReactionField;
 796             break;
 797
 798         case CoulombInteractionType::Switch:
 799         case CoulombInteractionType::Shift:
 800         case CoulombInteractionType::User:
 801         case CoulombInteractionType::PmeSwitch:
 802         case CoulombInteractionType::PmeUser:
 803         case CoulombInteractionType::PmeUserSwitch:
 804             forcerec->nbkernel_elec_interaction = NbkernelElecType::CubicSplineTable;
 805             break;
 806
 807         case CoulombInteractionType::Pme:
 808         case CoulombInteractionType::P3mAD:
 809         case CoulombInteractionType::Ewald:
 810             forcerec->nbkernel_elec_interaction = NbkernelElecType::Ewald;
 811             break;
 812
 813         default:
 814             gmx_fatal(FARGS,
 815                       "Unsupported electrostatic interaction: %s",
 816                       enumValueToString(interactionConst->eeltype));
 817     }
 818     forcerec->nbkernel_elec_modifier = interactionConst->coulomb_modifier;
 819
 820     /* Vdw: Translate from mdp settings to kernel format */
 821     switch (interactionConst->vdwtype)
 822     {
 823         case VanDerWaalsType::Cut:
 824             if (forcerec->haveBuckingham)
 825             {
 826                 forcerec->nbkernel_vdw_interaction = NbkernelVdwType::Buckingham;
 827             }
 828             else
 829             {
 830                 forcerec->nbkernel_vdw_interaction = NbkernelVdwType::LennardJones;
 831             }
 832             break;
 833         case VanDerWaalsType::Pme:
 834             forcerec->nbkernel_vdw_interaction = NbkernelVdwType::LJEwald;
 835             break;
 836
 837         case VanDerWaalsType::Switch:
 838         case VanDerWaalsType::Shift:
 839         case VanDerWaalsType::User:
 840             forcerec->nbkernel_vdw_interaction = NbkernelVdwType::CubicSplineTable;
 841             break;
 842
 843         default:
 844             gmx_fatal(FARGS, "Unsupported vdw interaction: %s", enumValueToString(interactionConst->vdwtype));
 845     }
 846     forcerec->nbkernel_vdw_modifier = interactionConst->vdw_modifier;
 847
 848     if (!gmx_within_tol(interactionConst->reppow, 12.0, 10 * GMX_DOUBLE_EPS))
 849     {
 850         gmx_fatal(FARGS, "Only LJ repulsion power 12 is supported");
 851     }
 852     /* Older tpr files can contain Coulomb user tables with the Verlet cutoff-scheme,
 853      * while mdrun does not (and never did) support this.
 854      */
 855     if (EEL_USER(forcerec->ic->eeltype))
 856     {
 857         gmx_fatal(FARGS,
 858                   "Electrostatics type %s is currently not supported",
 859                   enumValueToString(inputrec.coulombtype));
 860     }
 861
 862     /* 1-4 interaction electrostatics */
 863     forcerec->fudgeQQ = mtop.ffparams.fudgeQQ;
 864
 865     // Multiple time stepping
 866     forcerec->useMts = inputrec.useMts;
 867
 868     if (forcerec->useMts)
 869     {
 870         GMX_ASSERT(gmx::checkMtsRequirements(inputrec).empty(),
 871                    "All MTS requirements should be met here");
 872     }
 873
 874     const bool haveDirectVirialContributionsFast =
 875             forcerec->forceProviders->hasForceProvider() || gmx_mtop_ftype_count(mtop, F_POSRES) > 0
 876             || gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0 || inputrec.nwall > 0 || inputrec.bPull
 877             || inputrec.bRot || inputrec.bIMD;
 878     const bool haveDirectVirialContributionsSlow =
 879             EEL_FULL(interactionConst->eeltype) || EVDW_PME(interactionConst->vdwtype);
 880     for (int i = 0; i < (forcerec->useMts ? 2 : 1); i++)
 881     {
 882         bool haveDirectVirialContributions =
 883                 (((!forcerec->useMts || i == 0) && haveDirectVirialContributionsFast)
 884                  || ((!forcerec->useMts || i == 1) && haveDirectVirialContributionsSlow));
 885         forcerec->forceHelperBuffers.emplace_back(haveDirectVirialContributions);
 886     }
 887
 888     if (forcerec->shift_vec.empty())
 889     {
 890         forcerec->shift_vec.resize(gmx::c_numShiftVectors);
 891     }
 892
 893     GMX_ASSERT(forcerec->nbfp.empty(), "The nonbonded force parameters should not be set up yet.");
 894     forcerec->ntype = mtop.ffparams.atnr;
 895     forcerec->nbfp  = makeNonBondedParameterLists(
 896             mtop.ffparams.atnr, mtop.ffparams.iparams, forcerec->haveBuckingham);
 897     if (EVDW_PME(interactionConst->vdwtype))
 898     {
 899         forcerec->ljpme_c6grid = makeLJPmeC6GridCorrectionParameters(
 900                 mtop.ffparams.atnr, mtop.ffparams.iparams, forcerec->ljpme_combination_rule);
 901     }
 902
 903     /* Copy the energy group exclusions */
 904     forcerec->egp_flags = inputrec.opts.egp_flags;
 905
 906     /* Van der Waals stuff */
 907     if ((interactionConst->vdwtype != VanDerWaalsType::Cut)
 908         && (interactionConst->vdwtype != VanDerWaalsType::User) && !forcerec->haveBuckingham)
 909     {
 910         if (interactionConst->rvdw_switch >= interactionConst->rvdw)
 911         {
 912             gmx_fatal(FARGS,
 913                       "rvdw_switch (%f) must be < rvdw (%f)",
 914                       interactionConst->rvdw_switch,
 915                       interactionConst->rvdw);
 916         }
 917         if (fplog)
 918         {
 919             fprintf(fplog,
 920                     "Using %s Lennard-Jones, switch between %g and %g nm\n",
 921                     (interactionConst->eeltype == CoulombInteractionType::Switch) ? "switched" : "shifted",
 922                     interactionConst->rvdw_switch,
 923                     interactionConst->rvdw);
 924         }
 925     }
 926
 927     if (forcerec->haveBuckingham && EVDW_PME(interactionConst->vdwtype))
 928     {
 929         gmx_fatal(FARGS, "LJ PME not supported with Buckingham");
 930     }
 931
 932     if (forcerec->haveBuckingham
 933         && (interactionConst->vdwtype == VanDerWaalsType::Shift
 934             || interactionConst->vdwtype == VanDerWaalsType::Switch))
 935     {
 936         gmx_fatal(FARGS, "Switch/shift interaction not supported with Buckingham");
 937     }
 938
 939     if (forcerec->haveBuckingham)
 940     {
 941         gmx_fatal(FARGS, "The Verlet cutoff-scheme does not (yet) support Buckingham");
 942     }
 943
 944     if (inputrec.implicit_solvent)
 945     {
 946         gmx_fatal(FARGS, "Implict solvation is no longer supported.");
 947     }
 948
 949
 950     /* This code automatically gives table length tabext without cut-off's,
 951      * in that case grompp should already have checked that we do not need
 952      * normal tables and we only generate tables for 1-4 interactions.
 953      */
 954     real rtab = inputrec.rlist + inputrec.tabext;
 955
 956     /* We want to use unmodified tables for 1-4 coulombic
 957      * interactions, so we must in general have an extra set of
 958      * tables. */
 959     if (gmx_mtop_ftype_count(mtop, F_LJ14) > 0 || gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0
 960         || gmx_mtop_ftype_count(mtop, F_LJC_PAIRS_NB) > 0)
 961     {
 962         forcerec->pairsTable = make_tables(fplog, interactionConst, tabpfn, rtab, GMX_MAKETABLES_14ONLY);
 963     }
 964
 965     /* Wall stuff */
 966     forcerec->nwall = inputrec.nwall;
 967     if (inputrec.nwall && inputrec.wall_type == WallType::Table)
 968     {
 969         make_wall_tables(fplog, inputrec, tabfn, &mtop.groups, forcerec);
 970     }
 971
 972     forcerec->fcdata = std::make_unique<t_fcdata>();
 973
 974     if (!tabbfnm.empty())
 975     {
 976         t_fcdata& fcdata = *forcerec->fcdata;
 977         // Need to catch std::bad_alloc
 978         // TODO Don't need to catch this here, when merging with master branch
 979         try
 980         {
 981             // TODO move these tables into a separate struct and store reference in ListedForces
 982             fcdata.bondtab = make_bonded_tables(fplog, F_TABBONDS, F_TABBONDSNC, mtop, tabbfnm, "b");
 983             fcdata.angletab = make_bonded_tables(fplog, F_TABANGLES, -1, mtop, tabbfnm, "a");
 984             fcdata.dihtab   = make_bonded_tables(fplog, F_TABDIHS, -1, mtop, tabbfnm, "d");
 985         }
 986         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 987     }
 988     else
 989     {
 990         if (debug)
 991         {
 992             fprintf(debug,
 993                     "No fcdata or table file name passed, can not read table, can not do bonded "
 994                     "interactions\n");
 995         }
 996     }
 997
 998     /* Initialize the thread working data for bonded interactions */
 999     if (forcerec->useMts)
1000     {
1001         // Add one ListedForces object for each MTS level
1002         bool isFirstLevel = true;
1003         for (const auto& mtsLevel : inputrec.mtsLevels)
1004         {
1005             ListedForces::InteractionSelection interactionSelection;
1006             const auto&                        forceGroups = mtsLevel.forceGroups;
1007             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Pair)])
1008             {
1009                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Pairs));
1010             }
1011             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Dihedral)])
1012             {
1013                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Dihedrals));
1014             }
1015             if (forceGroups[static_cast<int>(gmx::MtsForceGroups::Angle)])
1016             {
1017                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Angles));
1018             }
1019             if (isFirstLevel)
1020             {
1021                 interactionSelection.set(static_cast<int>(ListedForces::InteractionGroup::Rest));
1022                 isFirstLevel = false;
1023             }
1024             forcerec->listedForces.emplace_back(
1025                     mtop.ffparams,
1026                     mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1027                     gmx_omp_nthreads_get(ModuleMultiThread::Bonded),
1028                     interactionSelection,
1029                     fplog);
1030         }
1031     }
1032     else
1033     {
1034         // Add one ListedForces object with all listed interactions
1035         forcerec->listedForces.emplace_back(
1036                 mtop.ffparams,
1037                 mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1038                 gmx_omp_nthreads_get(ModuleMultiThread::Bonded),
1039                 ListedForces::interactionSelectionAll(),
1040                 fplog);
1041     }
1042
1043     // QM/MM initialization if requested
1044     if (inputrec.bQMMM)
1045     {
1046         gmx_incons("QM/MM was requested, but is no longer available in GROMACS");
1047     }
1048
1049     /* Set all the static charge group info */
1050     forcerec->cginfo_mb = init_cginfo_mb(mtop, forcerec);
1051     if (!DOMAINDECOMP(commrec))
1052     {
1053         forcerec->cginfo = cginfo_expand(mtop.molblock.size(), forcerec->cginfo_mb);
1054     }
1055
1056     if (!DOMAINDECOMP(commrec))
1057     {
1058         forcerec_set_ranges(forcerec, mtop.natoms, mtop.natoms, mtop.natoms);
1059     }
1060
1061     forcerec->print_force = print_force;
1062
1063     forcerec->nthread_ewc = gmx_omp_nthreads_get(ModuleMultiThread::Bonded);
1064     forcerec->ewc_t.resize(forcerec->nthread_ewc);
1065
1066     if (inputrec.eDispCorr != DispersionCorrectionType::No)
1067     {
1068         forcerec->dispersionCorrection = std::make_unique<DispersionCorrection>(
1069                 mtop, inputrec, forcerec->haveBuckingham, forcerec->ntype, forcerec->nbfp, *forcerec->ic, tabfn);
1070         forcerec->dispersionCorrection->print(mdlog);
1071     }
1072
1073     if (fplog != nullptr)
1074     {
1075         /* Here we switch from using mdlog, which prints the newline before
1076          * the paragraph, to our old fprintf logging, which prints the newline
1077          * after the paragraph, so we should add a newline here.
1078          */
1079         fprintf(fplog, "\n");
1080     }
1081 }
1082
1083 t_forcerec::t_forcerec() = default;
1084
1085 t_forcerec::~t_forcerec() = default;