src/gromacs/gmxpreprocess/topio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "gromacs/utility/enumerationhelpers.h"
  41 #include "gromacs/utility/stringutil.h"
  42 #include "topio.h"
  43
  44 #include <cassert>
  45 #include <cctype>
  46 #include <cerrno>
  47 #include <cmath>
  48 #include <cstdio>
  49 #include <cstdlib>
  50 #include <cstring>
  51
  52 #include <algorithm>
  53 #include <memory>
  54 #include <unordered_set>
  55
  56 #include <sys/types.h>
  57
  58 #include "gromacs/fileio/gmxfio.h"
  59 #include "gromacs/fileio/warninp.h"
  60 #include "gromacs/gmxpreprocess/gmxcpp.h"
  61 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  62 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
  63 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
  64 #include "gromacs/gmxpreprocess/grompp_impl.h"
  65 #include "gromacs/gmxpreprocess/readir.h"
  66 #include "gromacs/gmxpreprocess/topdirs.h"
  67 #include "gromacs/gmxpreprocess/toppush.h"
  68 #include "gromacs/gmxpreprocess/topshake.h"
  69 #include "gromacs/gmxpreprocess/toputil.h"
  70 #include "gromacs/gmxpreprocess/vsite_parm.h"
  71 #include "gromacs/math/units.h"
  72 #include "gromacs/math/utilities.h"
  73 #include "gromacs/mdtypes/inputrec.h"
  74 #include "gromacs/mdtypes/md_enums.h"
  75 #include "gromacs/pbcutil/pbc.h"
  76 #include "gromacs/topology/block.h"
  77 #include "gromacs/topology/exclusionblocks.h"
  78 #include "gromacs/topology/ifunc.h"
  79 #include "gromacs/topology/symtab.h"
  80 #include "gromacs/topology/topology.h"
  81 #include "gromacs/utility/cstringutil.h"
  82 #include "gromacs/utility/fatalerror.h"
  83 #include "gromacs/utility/futil.h"
  84 #include "gromacs/utility/gmxassert.h"
  85 #include "gromacs/utility/logger.h"
  86 #include "gromacs/utility/pleasecite.h"
  87 #include "gromacs/utility/smalloc.h"
  88
  89 #define OPENDIR '['  /* starting sign for directive */
  90 #define CLOSEDIR ']' /* ending sign for directive   */
  91
  92 static void gen_pairs(const InteractionsOfType& nbs, InteractionsOfType* pairs, real fudge, CombinationRule comb)
  93 {
  94     real scaling;
  95     int  ntp = nbs.size();
  96     int  nnn = static_cast<int>(std::sqrt(static_cast<double>(ntp)));
  97     GMX_ASSERT(nnn * nnn == ntp,
  98                "Number of pairs of generated non-bonded parameters should be a perfect square");
  99     int nrfp  = NRFP(F_LJ);
 100     int nrfpA = interaction_function[F_LJ14].nrfpA;
 101     int nrfpB = interaction_function[F_LJ14].nrfpB;
 102
 103     if ((nrfp != nrfpA) || (nrfpA != nrfpB))
 104     {
 105         gmx_incons("Number of force parameters in gen_pairs wrong");
 106     }
 107
 108     fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
 109     pairs->interactionTypes.clear();
 110     int                             i = 0;
 111     std::array<int, 2>              atomNumbers;
 112     std::array<real, MAXFORCEPARAM> forceParam = { NOTSET };
 113     for (const auto& type : nbs.interactionTypes)
 114     {
 115         /* Copy type.atoms */
 116         atomNumbers = { i / nnn, i % nnn };
 117         /* Copy normal and FEP parameters and multiply by fudge factor */
 118         gmx::ArrayRef<const real> existingParam = type.forceParam();
 119         GMX_RELEASE_ASSERT(2 * nrfp <= MAXFORCEPARAM,
 120                            "Can't have more parameters than half of maximum p  arameter number");
 121         for (int j = 0; j < nrfp; j++)
 122         {
 123             /* If we are using sigma/epsilon values, only the epsilon values
 124              * should be scaled, but not sigma.
 125              * The sigma values have even indices 0,2, etc.
 126              */
 127             if ((comb == CombinationRule::Arithmetic || comb == CombinationRule::GeomSigEps)
 128                 && (j % 2 == 0))
 129             {
 130                 scaling = 1.0;
 131             }
 132             else
 133             {
 134                 scaling = fudge;
 135             }
 136
 137             forceParam[j]        = scaling * existingParam[j];
 138             forceParam[nrfp + j] = scaling * existingParam[j];
 139         }
 140         pairs->interactionTypes.emplace_back(InteractionOfType(atomNumbers, forceParam));
 141         i++;
 142     }
 143 }
 144
 145 double check_mol(const gmx_mtop_t* mtop, warninp* wi)
 146 {
 147     char   buf[256];
 148     int    i, ri;
 149     double q;
 150     real   m, mB;
 151
 152     /* Check mass and charge */
 153     q = 0.0;
 154
 155     for (const gmx_molblock_t& molb : mtop->molblock)
 156     {
 157         const t_atoms* atoms = &mtop->moltype[molb.type].atoms;
 158         for (i = 0; (i < atoms->nr); i++)
 159         {
 160             q += molb.nmol * atoms->atom[i].q;
 161             m               = atoms->atom[i].m;
 162             mB              = atoms->atom[i].mB;
 163             ParticleType pt = atoms->atom[i].ptype;
 164             /* If the particle is an atom or a nucleus it must have a mass,
 165              * else, if it is a shell, a vsite or a bondshell it can have mass zero
 166              */
 167             if (((m <= 0.0) || (mB <= 0.0)) && ((pt == ParticleType::Atom) || (pt == ParticleType::Nucleus)))
 168             {
 169                 ri = atoms->atom[i].resind;
 170                 sprintf(buf,
 171                         "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
 172                         *(atoms->atomname[i]),
 173                         *(atoms->resinfo[ri].name),
 174                         atoms->resinfo[ri].nr,
 175                         m,
 176                         mB);
 177                 warning_error(wi, buf);
 178             }
 179             else if (((m != 0) || (mB != 0)) && (pt == ParticleType::VSite))
 180             {
 181                 ri = atoms->atom[i].resind;
 182                 sprintf(buf,
 183                         "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state "
 184                         "B)\n"
 185                         "     Check your topology.\n",
 186                         *(atoms->atomname[i]),
 187                         *(atoms->resinfo[ri].name),
 188                         atoms->resinfo[ri].nr,
 189                         m,
 190                         mB);
 191                 warning_error(wi, buf);
 192                 /* The following statements make LINCS break! */
 193                 /* atoms->atom[i].m=0; */
 194             }
 195         }
 196     }
 197     return q;
 198 }
 199
 200 /*! \brief Returns the rounded charge of a molecule, when close to integer, otherwise returns the original charge.
 201  *
 202  * The results of this routine are only used for checking and for
 203  * printing warning messages. Thus we can assume that charges of molecules
 204  * should be integer. If the user wanted non-integer molecular charge,
 205  * an undesired warning is printed and the user should use grompp -maxwarn 1.
 206  *
 207  * \param qMol     The total, unrounded, charge of the molecule
 208  * \param sumAbsQ  The sum of absolute values of the charges, used for determining the tolerance for the rounding.
 209  */
 210 static double roundedMoleculeCharge(double qMol, double sumAbsQ)
 211 {
 212     /* We use a tolerance of 1e-6 for inaccuracies beyond the 6th decimal
 213      * of the charges for ascii float truncation in the topology files.
 214      * Although the summation here uses double precision, the charges
 215      * are read and stored in single precision when real=float. This can
 216      * lead to rounding errors of half the least significant bit.
 217      * Note that, unfortunately, we can not assume addition of random
 218      * rounding errors. It is not entirely unlikely that many charges
 219      * have a near half-bit rounding error with the same sign.
 220      */
 221     double tolAbs = 1e-6;
 222     double tol    = std::max(tolAbs, 0.5 * GMX_REAL_EPS * sumAbsQ);
 223     double qRound = std::round(qMol);
 224     if (std::abs(qMol - qRound) <= tol)
 225     {
 226         return qRound;
 227     }
 228     else
 229     {
 230         return qMol;
 231     }
 232 }
 233
 234 static void sum_q(const t_atoms* atoms, int numMols, double* qTotA, double* qTotB)
 235 {
 236     /* sum charge */
 237     double qmolA    = 0;
 238     double qmolB    = 0;
 239     double sumAbsQA = 0;
 240     double sumAbsQB = 0;
 241     for (int i = 0; i < atoms->nr; i++)
 242     {
 243         qmolA += atoms->atom[i].q;
 244         qmolB += atoms->atom[i].qB;
 245         sumAbsQA += std::abs(atoms->atom[i].q);
 246         sumAbsQB += std::abs(atoms->atom[i].qB);
 247     }
 248
 249     *qTotA += numMols * roundedMoleculeCharge(qmolA, sumAbsQA);
 250     *qTotB += numMols * roundedMoleculeCharge(qmolB, sumAbsQB);
 251 }
 252
 253 static void get_nbparm(char* nb_str, char* comb_str, VanDerWaalsPotential* nb, CombinationRule* comb, warninp* wi)
 254 {
 255     *nb = VanDerWaalsPotential::Count;
 256     for (auto i : gmx::EnumerationArray<VanDerWaalsPotential, bool>::keys())
 257     {
 258         if (gmx_strcasecmp(nb_str, enumValueToString(i)) == 0)
 259         {
 260             *nb = i;
 261         }
 262     }
 263     if (*nb == VanDerWaalsPotential::Count)
 264     {
 265         int integerValue = strtol(nb_str, nullptr, 10);
 266         if ((integerValue < 1) || (integerValue >= static_cast<int>(VanDerWaalsPotential::Count)))
 267         {
 268             std::string message =
 269                     gmx::formatString("Invalid nonbond function selector '%s' using %s",
 270                                       nb_str,
 271                                       enumValueToString(VanDerWaalsPotential::LJ));
 272             warning_error(wi, message);
 273             *nb = VanDerWaalsPotential::LJ;
 274         }
 275         else
 276         {
 277             *nb = static_cast<VanDerWaalsPotential>(integerValue);
 278         }
 279     }
 280     *comb = CombinationRule::Count;
 281     for (auto i : gmx::EnumerationArray<CombinationRule, bool>::keys())
 282     {
 283         if (gmx_strcasecmp(comb_str, enumValueToString(i)) == 0)
 284         {
 285             *comb = i;
 286         }
 287     }
 288     if (*comb == CombinationRule::Count)
 289     {
 290         int integerValue = strtol(comb_str, nullptr, 10);
 291         if ((integerValue < 1) || (integerValue >= static_cast<int>(CombinationRule::Count)))
 292         {
 293             std::string message =
 294                     gmx::formatString("Invalid combination rule selector '%s' using %s",
 295                                       comb_str,
 296                                       enumValueToString(CombinationRule::Geometric));
 297             warning_error(wi, message);
 298             *comb = CombinationRule::Geometric;
 299         }
 300         else
 301         {
 302             *comb = static_cast<CombinationRule>(integerValue);
 303         }
 304     }
 305 }
 306
 307 /*! \brief Parses define and include flags.
 308  *
 309  * Returns a vector of parsed include/define flags, with an extra nullptr entry at the back
 310  * for consumers that expect null-terminated char** structures.
 311  */
 312 static std::vector<char*> cpp_opts(const char* define, const char* include, warninp* wi)
 313 {
 314     int         n, len;
 315     const char* cppadds[2];
 316     const char* option[2] = { "-D", "-I" };
 317     const char* nopt[2]   = { "define", "include" };
 318     const char* ptr;
 319     const char* rptr;
 320     char*       buf;
 321     char        warn_buf[STRLEN];
 322
 323     cppadds[0] = define;
 324     cppadds[1] = include;
 325     std::vector<char*> cppOptions;
 326     for (n = 0; (n < 2); n++)
 327     {
 328         if (cppadds[n])
 329         {
 330             ptr = cppadds[n];
 331             while (*ptr != '\0')
 332             {
 333                 while ((*ptr != '\0') && isspace(*ptr))
 334                 {
 335                     ptr++;
 336                 }
 337                 rptr = ptr;
 338                 while ((*rptr != '\0') && !isspace(*rptr))
 339                 {
 340                     rptr++;
 341                 }
 342                 len = (rptr - ptr);
 343                 if (len > 2)
 344                 {
 345                     snew(buf, (len + 1));
 346                     strncpy(buf, ptr, len);
 347                     if (strstr(ptr, option[n]) != ptr)
 348                     {
 349                         set_warning_line(wi, "mdp file", -1);
 350                         sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
 351                         warning(wi, warn_buf);
 352                     }
 353                     else
 354                     {
 355                         cppOptions.emplace_back(gmx_strdup(buf));
 356                     }
 357                     sfree(buf);
 358                     ptr = rptr;
 359                 }
 360             }
 361         }
 362     }
 363     // Users of cppOptions expect a null last element.
 364     cppOptions.emplace_back(nullptr);
 365     return cppOptions;
 366 }
 367
 368
 369 static void make_atoms_sys(gmx::ArrayRef<const gmx_molblock_t>      molblock,
 370                            gmx::ArrayRef<const MoleculeInformation> molinfo,
 371                            t_atoms*                                 atoms)
 372 {
 373     atoms->nr   = 0;
 374     atoms->atom = nullptr;
 375
 376     for (const gmx_molblock_t& molb : molblock)
 377     {
 378         const t_atoms& mol_atoms = molinfo[molb.type].atoms;
 379
 380         srenew(atoms->atom, atoms->nr + molb.nmol * mol_atoms.nr);
 381
 382         for (int m = 0; m < molb.nmol; m++)
 383         {
 384             for (int a = 0; a < mol_atoms.nr; a++)
 385             {
 386                 atoms->atom[atoms->nr++] = mol_atoms.atom[a];
 387             }
 388         }
 389     }
 390 }
 391
 392
 393 static char** read_topol(const char*                           infile,
 394                          const char*                           outfile,
 395                          const char*                           define,
 396                          const char*                           include,
 397                          t_symtab*                             symtab,
 398                          PreprocessingAtomTypes*               atypes,
 399                          std::vector<MoleculeInformation>*     molinfo,
 400                          std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
 401                          gmx::ArrayRef<InteractionsOfType>     interactions,
 402                          CombinationRule*                      combination_rule,
 403                          double*                               reppow,
 404                          t_gromppopts*                         opts,
 405                          real*                                 fudgeQQ,
 406                          std::vector<gmx_molblock_t>*          molblock,
 407                          bool*                                 ffParametrizedWithHBondConstraints,
 408                          bool                                  bFEP,
 409                          bool                                  bZero,
 410                          bool                                  usingFullRangeElectrostatics,
 411                          warninp*                              wi,
 412                          const gmx::MDLogger&                  logger)
 413 {
 414     FILE*                out;
 415     int                  sl;
 416     char *               pline = nullptr, **title = nullptr;
 417     char                 line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
 418     char                 genpairs[32];
 419     char *               dirstr, *dummy2;
 420     int                  nrcopies, nscan, ncombs, ncopy;
 421     double               fLJ, fQQ, fPOW;
 422     MoleculeInformation* mi0 = nullptr;
 423     DirStack*            DS;
 424     Directive            d, newd;
 425     t_nbparam **         nbparam, **pair;
 426     real                 fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
 427     bool                 bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
 428     double               qt = 0, qBt = 0; /* total charge */
 429     int                  dcatt = -1, nmol_couple;
 430     /* File handling variables */
 431     int         status;
 432     bool        done;
 433     gmx_cpp_t   handle;
 434     char*       tmp_line = nullptr;
 435     char        warn_buf[STRLEN];
 436     const char* floating_point_arithmetic_tip =
 437             "Total charge should normally be an integer. See\n"
 438             "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
 439             "for discussion on how close it should be to an integer.\n";
 440     /* We need to open the output file before opening the input file,
 441      * because cpp_open_file can change the current working directory.
 442      */
 443     if (outfile)
 444     {
 445         out = gmx_fio_fopen(outfile, "w");
 446     }
 447     else
 448     {
 449         out = nullptr;
 450     }
 451
 452     /* open input file */
 453     auto cpp_opts_return = cpp_opts(define, include, wi);
 454     status               = cpp_open_file(infile, &handle, cpp_opts_return.data());
 455     if (status != 0)
 456     {
 457         gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
 458     }
 459
 460     /* some local variables */
 461     DS_Init(&DS);                   /* directive stack  */
 462     d       = Directive::d_invalid; /* first thing should be a directive */
 463     nbparam = nullptr;              /* The temporary non-bonded matrix */
 464     pair    = nullptr;              /* The temporary pair interaction matrix */
 465     std::vector<std::vector<gmx::ExclusionBlock>> exclusionBlocks;
 466     VanDerWaalsPotential                          nb_funct = VanDerWaalsPotential::LJ;
 467
 468     *reppow = 12.0; /* Default value for repulsion power     */
 469
 470     /* Init the number of CMAP torsion angles  and grid spacing */
 471     interactions[F_CMAP].cmakeGridSpacing = 0;
 472     interactions[F_CMAP].cmapAngles       = 0;
 473
 474     bWarn_copy_A_B = bFEP;
 475
 476     PreprocessingBondAtomType bondAtomType;
 477     /* parse the actual file */
 478     bReadDefaults = FALSE;
 479     bGenPairs     = FALSE;
 480     bReadMolType  = FALSE;
 481     nmol_couple   = 0;
 482
 483     do
 484     {
 485         status = cpp_read_line(&handle, STRLEN, line);
 486         done   = (status == eCPP_EOF);
 487         if (!done)
 488         {
 489             if (status != eCPP_OK)
 490             {
 491                 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
 492             }
 493             else if (out)
 494             {
 495                 fprintf(out, "%s\n", line);
 496             }
 497
 498             set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
 499
 500             pline = gmx_strdup(line);
 501
 502             /* Strip trailing '\' from pline, if it exists */
 503             sl = strlen(pline);
 504             if ((sl > 0) && (pline[sl - 1] == CONTINUE))
 505             {
 506                 pline[sl - 1] = ' ';
 507             }
 508
 509             /* build one long line from several fragments - necessary for CMAP */
 510             while (continuing(line))
 511             {
 512                 status = cpp_read_line(&handle, STRLEN, line);
 513                 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
 514
 515                 /* Since we depend on the '\' being present to continue to read, we copy line
 516                  * to a tmp string, strip the '\' from that string, and cat it to pline
 517                  */
 518                 tmp_line = gmx_strdup(line);
 519
 520                 sl = strlen(tmp_line);
 521                 if ((sl > 0) && (tmp_line[sl - 1] == CONTINUE))
 522                 {
 523                     tmp_line[sl - 1] = ' ';
 524                 }
 525
 526                 done = (status == eCPP_EOF);
 527                 if (!done)
 528                 {
 529                     if (status != eCPP_OK)
 530                     {
 531                         gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
 532                     }
 533                     else if (out)
 534                     {
 535                         fprintf(out, "%s\n", line);
 536                     }
 537                 }
 538
 539                 srenew(pline, strlen(pline) + strlen(tmp_line) + 1);
 540                 strcat(pline, tmp_line);
 541                 sfree(tmp_line);
 542             }
 543
 544             /* skip trailing and leading spaces and comment text */
 545             strip_comment(pline);
 546             trim(pline);
 547
 548             /* if there is something left... */
 549             if (static_cast<int>(strlen(pline)) > 0)
 550             {
 551                 if (pline[0] == OPENDIR)
 552                 {
 553                     /* A directive on this line: copy the directive
 554                      * without the brackets into dirstr, then
 555                      * skip spaces and tabs on either side of directive
 556                      */
 557                     dirstr = gmx_strdup((pline + 1));
 558                     if ((dummy2 = strchr(dirstr, CLOSEDIR)) != nullptr)
 559                     {
 560                         (*dummy2) = 0;
 561                     }
 562                     trim(dirstr);
 563
 564                     if ((newd = str2dir(dirstr)) == Directive::d_invalid)
 565                     {
 566                         sprintf(errbuf, "Invalid directive %s", dirstr);
 567                         warning_error(wi, errbuf);
 568                     }
 569                     else
 570                     {
 571                         /* Directive found */
 572                         if (DS_Check_Order(DS, newd))
 573                         {
 574                             DS_Push(&DS, newd);
 575                             d = newd;
 576                         }
 577                         else
 578                         {
 579                             /* we should print here which directives should have
 580                                been present, and which actually are */
 581                             gmx_fatal(FARGS,
 582                                       "%s\nInvalid order for directive %s",
 583                                       cpp_error(&handle, eCPP_SYNTAX),
 584                                       enumValueToString(newd));
 585                             /* d = Directive::d_invalid; */
 586                         }
 587
 588                         if (d == Directive::d_intermolecular_interactions)
 589                         {
 590                             if (*intermolecular_interactions == nullptr)
 591                             {
 592                                 /* We (mis)use the moleculetype processing
 593                                  * to process the intermolecular interactions
 594                                  * by making a "molecule" of the size of the system.
 595                                  */
 596                                 *intermolecular_interactions = std::make_unique<MoleculeInformation>();
 597                                 mi0                          = intermolecular_interactions->get();
 598                                 mi0->initMolInfo();
 599                                 make_atoms_sys(*molblock, *molinfo, &mi0->atoms);
 600                             }
 601                         }
 602                     }
 603                     sfree(dirstr);
 604                 }
 605                 else if (d != Directive::d_invalid)
 606                 {
 607                     /* Not a directive, just a plain string
 608                      * use a gigantic switch to decode,
 609                      * if there is a valid directive!
 610                      */
 611                     switch (d)
 612                     {
 613                         case Directive::d_defaults:
 614                             if (bReadDefaults)
 615                             {
 616                                 gmx_fatal(FARGS,
 617                                           "%s\nFound a second defaults directive.\n",
 618                                           cpp_error(&handle, eCPP_SYNTAX));
 619                             }
 620                             bReadDefaults = TRUE;
 621                             nscan         = sscanf(
 622                                     pline, "%s%s%s%lf%lf%lf", nb_str, comb_str, genpairs, &fLJ, &fQQ, &fPOW);
 623                             if (nscan < 2)
 624                             {
 625                                 too_few(wi);
 626                             }
 627                             else
 628                             {
 629                                 bGenPairs = FALSE;
 630                                 fudgeLJ   = 1.0;
 631                                 *fudgeQQ  = 1.0;
 632
 633                                 get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
 634                                 if (nscan >= 3)
 635                                 {
 636                                     bGenPairs = (gmx::equalCaseInsensitive(genpairs, "Y", 1));
 637                                     if (nb_funct != VanDerWaalsPotential::LJ && bGenPairs)
 638                                     {
 639                                         gmx_fatal(FARGS,
 640                                                   "Generating pair parameters is only supported "
 641                                                   "with LJ non-bonded interactions");
 642                                     }
 643                                 }
 644                                 if (nscan >= 4)
 645                                 {
 646                                     fudgeLJ = fLJ;
 647                                 }
 648                                 if (nscan >= 5)
 649                                 {
 650                                     *fudgeQQ = fQQ;
 651                                 }
 652                                 if (nscan >= 6)
 653                                 {
 654                                     *reppow = fPOW;
 655                                 }
 656                             }
 657                             nb_funct = static_cast<VanDerWaalsPotential>(ifunc_index(
 658                                     Directive::d_nonbond_params, static_cast<int>(nb_funct)));
 659
 660                             break;
 661                         case Directive::d_atomtypes:
 662                             push_at(symtab,
 663                                     atypes,
 664                                     &bondAtomType,
 665                                     pline,
 666                                     static_cast<int>(nb_funct),
 667                                     &nbparam,
 668                                     bGenPairs ? &pair : nullptr,
 669                                     wi);
 670                             break;
 671
 672                         case Directive::d_bondtypes: // Intended to fall through
 673                         case Directive::d_constrainttypes:
 674                             push_bt(d, interactions, 2, nullptr, &bondAtomType, pline, wi);
 675                             break;
 676                         case Directive::d_pairtypes:
 677                             if (bGenPairs)
 678                             {
 679                                 push_nbt(d, pair, atypes, pline, F_LJ14, wi);
 680                             }
 681                             else
 682                             {
 683                                 push_bt(d, interactions, 2, atypes, nullptr, pline, wi);
 684                             }
 685                             break;
 686                         case Directive::d_angletypes:
 687                             push_bt(d, interactions, 3, nullptr, &bondAtomType, pline, wi);
 688                             break;
 689                         case Directive::d_dihedraltypes:
 690                             /* Special routine that can read both 2 and 4 atom dihedral definitions. */
 691                             push_dihedraltype(d, interactions, &bondAtomType, pline, wi);
 692                             break;
 693
 694                         case Directive::d_nonbond_params:
 695                             push_nbt(d, nbparam, atypes, pline, static_cast<int>(nb_funct), wi);
 696                             break;
 697
 698                         case Directive::d_implicit_genborn_params: // NOLINT bugprone-branch-clone
 699                             // Skip this line, so old topologies with
 700                             // GB parameters can be read.
 701                             break;
 702
 703                         case Directive::d_implicit_surface_params:
 704                             // Skip this line, so that any topologies
 705                             // with surface parameters can be read
 706                             // (even though these were never formally
 707                             // supported).
 708                             break;
 709
 710                         case Directive::d_cmaptypes:
 711                             push_cmaptype(d, interactions, 5, atypes, &bondAtomType, pline, wi);
 712                             break;
 713
 714                         case Directive::d_moleculetype:
 715                         {
 716                             if (!bReadMolType)
 717                             {
 718                                 int ntype;
 719                                 if (opts->couple_moltype != nullptr
 720                                     && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam0 == ecouplamQ
 721                                         || opts->couple_lam1 == ecouplamNONE
 722                                         || opts->couple_lam1 == ecouplamQ))
 723                                 {
 724                                     dcatt = add_atomtype_decoupled(
 725                                             symtab, atypes, &nbparam, bGenPairs ? &pair : nullptr);
 726                                 }
 727                                 ntype  = atypes->size();
 728                                 ncombs = (ntype * (ntype + 1)) / 2;
 729                                 generate_nbparams(*combination_rule,
 730                                                   static_cast<int>(nb_funct),
 731                                                   &(interactions[static_cast<int>(nb_funct)]),
 732                                                   atypes,
 733                                                   wi);
 734                                 ncopy = copy_nbparams(nbparam,
 735                                                       static_cast<int>(nb_funct),
 736                                                       &(interactions[static_cast<int>(nb_funct)]),
 737                                                       ntype);
 738                                 GMX_LOG(logger.info)
 739                                         .asParagraph()
 740                                         .appendTextFormatted(
 741                                                 "Generated %d of the %d non-bonded parameter "
 742                                                 "combinations",
 743                                                 ncombs - ncopy,
 744                                                 ncombs);
 745                                 free_nbparam(nbparam, ntype);
 746                                 if (bGenPairs)
 747                                 {
 748                                     gen_pairs((interactions[static_cast<int>(nb_funct)]),
 749                                               &(interactions[F_LJ14]),
 750                                               fudgeLJ,
 751                                               *combination_rule);
 752                                     ncopy = copy_nbparams(
 753                                             pair, static_cast<int>(nb_funct), &(interactions[F_LJ14]), ntype);
 754                                     GMX_LOG(logger.info)
 755                                             .asParagraph()
 756                                             .appendTextFormatted(
 757                                                     "Generated %d of the %d 1-4 parameter "
 758                                                     "combinations",
 759                                                     ncombs - ncopy,
 760                                                     ncombs);
 761                                     free_nbparam(pair, ntype);
 762                                 }
 763                                 /* Copy GBSA parameters to atomtype array? */
 764
 765                                 bReadMolType = TRUE;
 766                             }
 767
 768                             push_molt(symtab, molinfo, pline, wi);
 769                             exclusionBlocks.emplace_back();
 770                             mi0                    = &molinfo->back();
 771                             mi0->atoms.haveMass    = TRUE;
 772                             mi0->atoms.haveCharge  = TRUE;
 773                             mi0->atoms.haveType    = TRUE;
 774                             mi0->atoms.haveBState  = TRUE;
 775                             mi0->atoms.havePdbInfo = FALSE;
 776                             break;
 777                         }
 778                         case Directive::d_atoms:
 779                             push_atom(symtab, &(mi0->atoms), atypes, pline, wi);
 780                             break;
 781
 782                         case Directive::d_pairs:
 783                             GMX_RELEASE_ASSERT(
 784                                     mi0,
 785                                     "Need to have a valid MoleculeInformation object to work on");
 786                             push_bond(d,
 787                                       interactions,
 788                                       mi0->interactions,
 789                                       &(mi0->atoms),
 790                                       atypes,
 791                                       pline,
 792                                       FALSE,
 793                                       bGenPairs,
 794                                       *fudgeQQ,
 795                                       bZero,
 796                                       &bWarn_copy_A_B,
 797                                       wi);
 798                             break;
 799                         case Directive::d_pairs_nb:
 800                             GMX_RELEASE_ASSERT(
 801                                     mi0,
 802                                     "Need to have a valid MoleculeInformation object to work on");
 803                             push_bond(d,
 804                                       interactions,
 805                                       mi0->interactions,
 806                                       &(mi0->atoms),
 807                                       atypes,
 808                                       pline,
 809                                       FALSE,
 810                                       FALSE,
 811                                       1.0,
 812                                       bZero,
 813                                       &bWarn_copy_A_B,
 814                                       wi);
 815                             break;
 816
 817                         case Directive::d_vsites1:
 818                         case Directive::d_vsites2:
 819                         case Directive::d_vsites3:
 820                         case Directive::d_vsites4:
 821                         case Directive::d_bonds:
 822                         case Directive::d_angles:
 823                         case Directive::d_constraints:
 824                         case Directive::d_settles:
 825                         case Directive::d_position_restraints:
 826                         case Directive::d_angle_restraints:
 827                         case Directive::d_angle_restraints_z:
 828                         case Directive::d_distance_restraints:
 829                         case Directive::d_orientation_restraints:
 830                         case Directive::d_dihedral_restraints:
 831                         case Directive::d_dihedrals:
 832                         case Directive::d_polarization:
 833                         case Directive::d_water_polarization:
 834                         case Directive::d_thole_polarization:
 835                             GMX_RELEASE_ASSERT(
 836                                     mi0,
 837                                     "Need to have a valid MoleculeInformation object to work on");
 838                             push_bond(d,
 839                                       interactions,
 840                                       mi0->interactions,
 841                                       &(mi0->atoms),
 842                                       atypes,
 843                                       pline,
 844                                       TRUE,
 845                                       bGenPairs,
 846                                       *fudgeQQ,
 847                                       bZero,
 848                                       &bWarn_copy_A_B,
 849                                       wi);
 850                             break;
 851                         case Directive::d_cmap:
 852                             GMX_RELEASE_ASSERT(
 853                                     mi0,
 854                                     "Need to have a valid MoleculeInformation object to work on");
 855                             push_cmap(d, interactions, mi0->interactions, &(mi0->atoms), atypes, pline, wi);
 856                             break;
 857
 858                         case Directive::d_vsitesn:
 859                             GMX_RELEASE_ASSERT(
 860                                     mi0,
 861                                     "Need to have a valid MoleculeInformation object to work on");
 862                             push_vsitesn(d, mi0->interactions, &(mi0->atoms), pline, wi);
 863                             break;
 864                         case Directive::d_exclusions:
 865                             GMX_ASSERT(!exclusionBlocks.empty(),
 866                                        "exclusionBlocks must always be allocated so exclusions can "
 867                                        "be processed");
 868                             if (exclusionBlocks.back().empty())
 869                             {
 870                                 GMX_RELEASE_ASSERT(mi0,
 871                                                    "Need to have a valid MoleculeInformation "
 872                                                    "object to work on");
 873                                 exclusionBlocks.back().resize(mi0->atoms.nr);
 874                             }
 875                             push_excl(pline, exclusionBlocks.back(), wi);
 876                             break;
 877                         case Directive::d_system:
 878                             trim(pline);
 879                             title = put_symtab(symtab, pline);
 880                             break;
 881                         case Directive::d_molecules:
 882                         {
 883                             int  whichmol;
 884                             bool bCouple;
 885
 886                             push_mol(*molinfo, pline, &whichmol, &nrcopies, wi);
 887                             mi0 = &((*molinfo)[whichmol]);
 888                             molblock->resize(molblock->size() + 1);
 889                             molblock->back().type = whichmol;
 890                             molblock->back().nmol = nrcopies;
 891
 892                             bCouple = (opts->couple_moltype != nullptr
 893                                        && (gmx_strcasecmp("system", opts->couple_moltype) == 0
 894                                            || strcmp(*(mi0->name), opts->couple_moltype) == 0));
 895                             if (bCouple)
 896                             {
 897                                 nmol_couple += nrcopies;
 898                             }
 899
 900                             if (mi0->atoms.nr == 0)
 901                             {
 902                                 gmx_fatal(FARGS, "Molecule type '%s' contains no atoms", *mi0->name);
 903                             }
 904                             GMX_LOG(logger.info)
 905                                     .asParagraph()
 906                                     .appendTextFormatted(
 907                                             "Excluding %d bonded neighbours molecule type '%s'",
 908                                             mi0->nrexcl,
 909                                             *mi0->name);
 910                             sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
 911                             if (!mi0->bProcessed)
 912                             {
 913                                 generate_excl(mi0->nrexcl, mi0->atoms.nr, mi0->interactions, &(mi0->excls));
 914                                 gmx::mergeExclusions(&(mi0->excls), exclusionBlocks[whichmol]);
 915                                 make_shake(mi0->interactions, &mi0->atoms, opts->nshake, logger);
 916
 917                                 if (bCouple)
 918                                 {
 919                                     convert_moltype_couple(mi0,
 920                                                            dcatt,
 921                                                            *fudgeQQ,
 922                                                            opts->couple_lam0,
 923                                                            opts->couple_lam1,
 924                                                            opts->bCoupleIntra,
 925                                                            static_cast<int>(nb_funct),
 926                                                            &(interactions[static_cast<int>(nb_funct)]),
 927                                                            wi);
 928                                 }
 929                                 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
 930                                 mi0->bProcessed = TRUE;
 931                             }
 932                             break;
 933                         }
 934                         default:
 935                             GMX_LOG(logger.warning)
 936                                     .asParagraph()
 937                                     .appendTextFormatted("case: %d", static_cast<int>(d));
 938                             gmx_incons("unknown directive");
 939                     }
 940                 }
 941             }
 942             sfree(pline);
 943             pline = nullptr;
 944         }
 945     } while (!done);
 946
 947     // Check that all strings defined with -D were used when processing topology
 948     std::string unusedDefineWarning = checkAndWarnForUnusedDefines(*handle);
 949     if (!unusedDefineWarning.empty())
 950     {
 951         warning(wi, unusedDefineWarning);
 952     }
 953
 954     for (char* element : cpp_opts_return)
 955     {
 956         sfree(element);
 957     }
 958
 959     if (out)
 960     {
 961         gmx_fio_fclose(out);
 962     }
 963
 964     /* List of GROMACS define names for force fields that have been
 965      * parametrized using constraints involving hydrogens only.
 966      *
 967      * We should avoid hardcoded names, but this is hopefully only
 968      * needed temparorily for discouraging use of constraints=all-bonds.
 969      */
 970     const std::array<std::string, 3> ffDefines = { "_FF_AMBER", "_FF_CHARMM", "_FF_OPLSAA" };
 971     *ffParametrizedWithHBondConstraints        = false;
 972     for (const std::string& ffDefine : ffDefines)
 973     {
 974         if (cpp_find_define(&handle, ffDefine))
 975         {
 976             *ffParametrizedWithHBondConstraints = true;
 977         }
 978     }
 979
 980     if (cpp_find_define(&handle, "_FF_GROMOS96") != nullptr)
 981     {
 982         warning(wi,
 983                 "The GROMOS force fields have been parametrized with a physically incorrect "
 984                 "multiple-time-stepping scheme for a twin-range cut-off. When used with "
 985                 "a single-range cut-off (or a correct Trotter multiple-time-stepping scheme), "
 986                 "physical properties, such as the density, might differ from the intended values. "
 987                 "Since there are researchers actively working on validating GROMOS with modern "
 988                 "integrators we have not yet removed the GROMOS force fields, but you should be "
 989                 "aware of these issues and check if molecules in your system are affected before "
 990                 "proceeding. "
 991                 "Further information is available at https://redmine.gromacs.org/issues/2884 , "
 992                 "and a longer explanation of our decision to remove physically incorrect "
 993                 "algorithms "
 994                 "can be found at https://doi.org/10.26434/chemrxiv.11474583.v1 .");
 995     }
 996     // TODO: Update URL for Issue #2884 in conjunction with updating grompp.warn in regressiontests.
 997
 998     cpp_done(handle);
 999
1000     if (opts->couple_moltype)
1001     {
1002         if (nmol_couple == 0)
1003         {
1004             gmx_fatal(FARGS, "Did not find any molecules of type '%s' for coupling", opts->couple_moltype);
1005         }
1006         GMX_LOG(logger.info)
1007                 .asParagraph()
1008                 .appendTextFormatted(
1009                         "Coupling %d copies of molecule type '%s'", nmol_couple, opts->couple_moltype);
1010     }
1011
1012     /* this is not very clean, but fixes core dump on empty system name */
1013     if (!title)
1014     {
1015         title = put_symtab(symtab, "");
1016     }
1017
1018     if (fabs(qt) > 1e-4)
1019     {
1020         sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
1021         warning_note(wi, warn_buf);
1022     }
1023     if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
1024     {
1025         sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt, floating_point_arithmetic_tip);
1026         warning_note(wi, warn_buf);
1027     }
1028     if (usingFullRangeElectrostatics && (fabs(qt) > 1e-4 || fabs(qBt) > 1e-4))
1029     {
1030         warning(wi,
1031                 "You are using Ewald electrostatics in a system with net charge. This can lead to "
1032                 "severe artifacts, such as ions moving into regions with low dielectric, due to "
1033                 "the uniform background charge. We suggest to neutralize your system with counter "
1034                 "ions, possibly in combination with a physiological salt concentration.");
1035         please_cite(stdout, "Hub2014a");
1036     }
1037
1038     DS_Done(&DS);
1039
1040     if (*intermolecular_interactions != nullptr)
1041     {
1042         sfree(intermolecular_interactions->get()->atoms.atom);
1043     }
1044
1045     return title;
1046 }
1047
1048 char** do_top(bool                                  bVerbose,
1049               const char*                           topfile,
1050               const char*                           topppfile,
1051               t_gromppopts*                         opts,
1052               bool                                  bZero,
1053               t_symtab*                             symtab,
1054               gmx::ArrayRef<InteractionsOfType>     interactions,
1055               CombinationRule*                      combination_rule,
1056               double*                               repulsion_power,
1057               real*                                 fudgeQQ,
1058               PreprocessingAtomTypes*               atypes,
1059               std::vector<MoleculeInformation>*     molinfo,
1060               std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
1061               const t_inputrec*                     ir,
1062               std::vector<gmx_molblock_t>*          molblock,
1063               bool*                                 ffParametrizedWithHBondConstraints,
1064               warninp*                              wi,
1065               const gmx::MDLogger&                  logger)
1066 {
1067     /* Tmpfile might contain a long path */
1068     const char* tmpfile;
1069     char**      title;
1070
1071     if (topppfile)
1072     {
1073         tmpfile = topppfile;
1074     }
1075     else
1076     {
1077         tmpfile = nullptr;
1078     }
1079
1080     if (bVerbose)
1081     {
1082         GMX_LOG(logger.info).asParagraph().appendTextFormatted("processing topology...");
1083     }
1084     title = read_topol(topfile,
1085                        tmpfile,
1086                        opts->define,
1087                        opts->include,
1088                        symtab,
1089                        atypes,
1090                        molinfo,
1091                        intermolecular_interactions,
1092                        interactions,
1093                        combination_rule,
1094                        repulsion_power,
1095                        opts,
1096                        fudgeQQ,
1097                        molblock,
1098                        ffParametrizedWithHBondConstraints,
1099                        ir->efep != FreeEnergyPerturbationType::No,
1100                        bZero,
1101                        EEL_FULL(ir->coulombtype),
1102                        wi,
1103                        logger);
1104
1105     if ((*combination_rule != CombinationRule::Geometric) && (ir->vdwtype == VanDerWaalsType::User))
1106     {
1107         warning(wi,
1108                 "Using sigma/epsilon based combination rules with"
1109                 " user supplied potential function may produce unwanted"
1110                 " results");
1111     }
1112
1113     return title;
1114 }
1115
1116 /*! \brief
1117  * Exclude molecular interactions for QM atoms handled by MiMic.
1118  *
1119  * Update the exclusion lists to include all QM atoms of this molecule,
1120  * replace bonds between QM atoms with CONNBOND and
1121  * set charges of QM atoms to 0.
1122  *
1123  * \param[in,out] molt molecule type with QM atoms
1124  * \param[in] grpnr group informatio
1125  * \param[in,out] ir input record
1126  * \param[in] logger Handle to logging interface.
1127  */
1128 static void generate_qmexcl_moltype(gmx_moltype_t*       molt,
1129                                     const unsigned char* grpnr,
1130                                     t_inputrec*          ir,
1131                                     const gmx::MDLogger& logger)
1132 {
1133     /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1134
1135     /* generates the exclusions between the individual QM atoms, as
1136      * these interactions should be handled by the QM subroutines and
1137      * not by the gromacs routines
1138      */
1139     int   qm_max = 0, qm_nr = 0, link_nr = 0;
1140     int * qm_arr = nullptr, *link_arr = nullptr;
1141     bool* bQMMM;
1142
1143     /* First we search and select the QM atoms in an qm_arr array that
1144      * we use to create the exclusions.
1145      *
1146      * we take the possibility into account that a user has defined more
1147      * than one QM group:
1148      *
1149      * for that we also need to do this an ugly work-about just in case
1150      * the QM group contains the entire system...
1151      */
1152
1153     /* we first search for all the QM atoms and put them in an array
1154      */
1155     for (int j = 0; j < ir->opts.ngQM; j++)
1156     {
1157         for (int i = 0; i < molt->atoms.nr; i++)
1158         {
1159             if (qm_nr >= qm_max)
1160             {
1161                 qm_max += 100;
1162                 srenew(qm_arr, qm_max);
1163             }
1164             if ((grpnr ? grpnr[i] : 0) == j)
1165             {
1166                 qm_arr[qm_nr++]        = i;
1167                 molt->atoms.atom[i].q  = 0.0;
1168                 molt->atoms.atom[i].qB = 0.0;
1169             }
1170         }
1171     }
1172     /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1173      * QM/not QM. We first set all elements to false. Afterwards we use
1174      * the qm_arr to change the elements corresponding to the QM atoms
1175      * to TRUE.
1176      */
1177     snew(bQMMM, molt->atoms.nr);
1178     for (int i = 0; i < molt->atoms.nr; i++)
1179     {
1180         bQMMM[i] = FALSE;
1181     }
1182     for (int i = 0; i < qm_nr; i++)
1183     {
1184         bQMMM[qm_arr[i]] = TRUE;
1185     }
1186
1187     /* We remove all bonded interactions (i.e. bonds,
1188      * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1189      * are removed is as follows: if the interaction invloves 2 atoms,
1190      * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1191      * it is removed if at least two of the atoms are QM atoms, if the
1192      * interaction involves 4 atoms, it is removed if there are at least
1193      * 2 QM atoms.  Since this routine is called once before any forces
1194      * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1195      * be rewritten at this poitn without any problem. 25-9-2002 */
1196
1197     /* first check whether we already have CONNBONDS.
1198      * Note that if we don't, we don't add a param entry and set ftype=0,
1199      * which is ok, since CONNBONDS does not use parameters.
1200      */
1201     int ftype_connbond = 0;
1202     int ind_connbond   = 0;
1203     if (!molt->ilist[F_CONNBONDS].empty())
1204     {
1205         GMX_LOG(logger.info)
1206                 .asParagraph()
1207                 .appendTextFormatted("nr. of CONNBONDS present already: %d",
1208                                      molt->ilist[F_CONNBONDS].size() / 3);
1209         ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0];
1210         ind_connbond   = molt->ilist[F_CONNBONDS].size();
1211     }
1212     /* now we delete all bonded interactions, except the ones describing
1213      * a chemical bond. These are converted to CONNBONDS
1214      */
1215     for (int ftype = 0; ftype < F_NRE; ftype++)
1216     {
1217         if (!(interaction_function[ftype].flags & IF_BOND) || ftype == F_CONNBONDS)
1218         {
1219             continue;
1220         }
1221         int nratoms = interaction_function[ftype].nratoms;
1222         int j       = 0;
1223         while (j < molt->ilist[ftype].size())
1224         {
1225             bool bexcl;
1226
1227             if (nratoms == 2)
1228             {
1229                 /* Remove an interaction between two atoms when both are
1230                  * in the QM region. Note that we don't have to worry about
1231                  * link atoms here, as they won't have 2-atom interactions.
1232                  */
1233                 int a1 = molt->ilist[ftype].iatoms[1 + j + 0];
1234                 int a2 = molt->ilist[ftype].iatoms[1 + j + 1];
1235                 bexcl  = (bQMMM[a1] && bQMMM[a2]);
1236                 /* A chemical bond between two QM atoms will be copied to
1237                  * the F_CONNBONDS list, for reasons mentioned above.
1238                  */
1239                 if (bexcl && IS_CHEMBOND(ftype))
1240                 {
1241                     InteractionList& ilist = molt->ilist[F_CONNBONDS];
1242                     ilist.iatoms.resize(ind_connbond + 3);
1243                     ilist.iatoms[ind_connbond++] = ftype_connbond;
1244                     ilist.iatoms[ind_connbond++] = a1;
1245                     ilist.iatoms[ind_connbond++] = a2;
1246                 }
1247             }
1248             else
1249             {
1250                 /* MM interactions have to be excluded if they are included
1251                  * in the QM already. Because we use a link atom (H atom)
1252                  * when the QM/MM boundary runs through a chemical bond, this
1253                  * means that as long as one atom is MM, we still exclude,
1254                  * as the interaction is included in the QM via:
1255                  * QMatom1-QMatom2-QMatom-3-Linkatom.
1256                  */
1257                 int numQmAtoms = 0;
1258                 for (int jj = j + 1; jj < j + 1 + nratoms; jj++)
1259                 {
1260                     if (bQMMM[molt->ilist[ftype].iatoms[jj]])
1261                     {
1262                         numQmAtoms++;
1263                     }
1264                 }
1265
1266                 /* MiMiC treats link atoms as quantum atoms - therefore
1267                  * we do not need do additional exclusions here */
1268                 bexcl = numQmAtoms == nratoms;
1269
1270                 if (bexcl && ftype == F_SETTLE)
1271                 {
1272                     gmx_fatal(FARGS,
1273                               "Can not apply QM to molecules with SETTLE, replace the moleculetype "
1274                               "using QM and SETTLE by one without SETTLE");
1275                 }
1276             }
1277             if (bexcl)
1278             {
1279                 /* since the interaction involves QM atoms, these should be
1280                  * removed from the MM ilist
1281                  */
1282                 InteractionList& ilist = molt->ilist[ftype];
1283                 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1284                 {
1285                     ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1286                 }
1287                 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1288             }
1289             else
1290             {
1291                 j += nratoms + 1; /* the +1 is for the functype */
1292             }
1293         }
1294     }
1295     /* Now, we search for atoms bonded to a QM atom because we also want
1296      * to exclude their nonbonded interactions with the QM atoms. The
1297      * reason for this is that this interaction is accounted for in the
1298      * linkatoms interaction with the QMatoms and would be counted
1299      * twice.  */
1300
1301     /* creating the exclusion block for the QM atoms. Each QM atom has
1302      * as excluded elements all the other QMatoms (and itself).
1303      */
1304     t_blocka qmexcl;
1305     qmexcl.nr  = molt->atoms.nr;
1306     qmexcl.nra = qm_nr * (qm_nr + link_nr) + link_nr * qm_nr;
1307     snew(qmexcl.index, qmexcl.nr + 1);
1308     snew(qmexcl.a, qmexcl.nra);
1309     int j = 0;
1310     for (int i = 0; i < qmexcl.nr; i++)
1311     {
1312         qmexcl.index[i] = j;
1313         if (bQMMM[i])
1314         {
1315             for (int k = 0; k < qm_nr; k++)
1316             {
1317                 qmexcl.a[k + j] = qm_arr[k];
1318             }
1319             for (int k = 0; k < link_nr; k++)
1320             {
1321                 qmexcl.a[qm_nr + k + j] = link_arr[k];
1322             }
1323             j += (qm_nr + link_nr);
1324         }
1325     }
1326     qmexcl.index[qmexcl.nr] = j;
1327
1328     /* and merging with the exclusions already present in sys.
1329      */
1330
1331     std::vector<gmx::ExclusionBlock> qmexcl2(molt->atoms.nr);
1332     gmx::blockaToExclusionBlocks(&qmexcl, qmexcl2);
1333     gmx::mergeExclusions(&(molt->excls), qmexcl2);
1334
1335     /* Finally, we also need to get rid of the pair interactions of the
1336      * classical atom bonded to the boundary QM atoms with the QMatoms,
1337      * as this interaction is already accounted for by the QM, so also
1338      * here we run the risk of double counting! We proceed in a similar
1339      * way as we did above for the other bonded interactions: */
1340     for (int i = F_LJ14; i < F_COUL14; i++)
1341     {
1342         int nratoms = interaction_function[i].nratoms;
1343         int j       = 0;
1344         while (j < molt->ilist[i].size())
1345         {
1346             int  a1    = molt->ilist[i].iatoms[j + 1];
1347             int  a2    = molt->ilist[i].iatoms[j + 2];
1348             bool bexcl = (bQMMM[a1] && bQMMM[a2]);
1349             if (bexcl)
1350             {
1351                 /* since the interaction involves QM atoms, these should be
1352                  * removed from the MM ilist
1353                  */
1354                 InteractionList& ilist = molt->ilist[i];
1355                 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1356                 {
1357                     ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1358                 }
1359                 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1360             }
1361             else
1362             {
1363                 j += nratoms + 1; /* the +1 is for the functype */
1364             }
1365         }
1366     }
1367
1368     free(qm_arr);
1369     free(bQMMM);
1370     free(link_arr);
1371 } /* generate_qmexcl */
1372
1373 void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, const gmx::MDLogger& logger)
1374 {
1375     /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1376      */
1377
1378     unsigned char*  grpnr;
1379     int             mol, nat_mol, nr_mol_with_qm_atoms = 0;
1380     gmx_molblock_t* molb;
1381     bool            bQMMM;
1382     int             index_offset = 0;
1383     int             qm_nr        = 0;
1384
1385     grpnr = sys->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].data();
1386
1387     for (size_t mb = 0; mb < sys->molblock.size(); mb++)
1388     {
1389         molb    = &sys->molblock[mb];
1390         nat_mol = sys->moltype[molb->type].atoms.nr;
1391         for (mol = 0; mol < molb->nmol; mol++)
1392         {
1393             bQMMM = FALSE;
1394             for (int i = 0; i < nat_mol; i++)
1395             {
1396                 if ((grpnr ? grpnr[i] : 0) < (ir->opts.ngQM))
1397                 {
1398                     bQMMM = TRUE;
1399                     qm_nr++;
1400                 }
1401             }
1402
1403             if (bQMMM)
1404             {
1405                 nr_mol_with_qm_atoms++;
1406                 if (molb->nmol > 1)
1407                 {
1408                     /* We need to split this molblock */
1409                     if (mol > 0)
1410                     {
1411                         /* Split the molblock at this molecule */
1412                         auto pos = sys->molblock.begin() + mb + 1;
1413                         sys->molblock.insert(pos, sys->molblock[mb]);
1414                         sys->molblock[mb].nmol = mol;
1415                         sys->molblock[mb + 1].nmol -= mol;
1416                         mb++;
1417                         molb = &sys->molblock[mb];
1418                     }
1419                     if (molb->nmol > 1)
1420                     {
1421                         /* Split the molblock after this molecule */
1422                         auto pos = sys->molblock.begin() + mb + 1;
1423                         sys->molblock.insert(pos, sys->molblock[mb]);
1424                         molb                   = &sys->molblock[mb];
1425                         sys->molblock[mb].nmol = 1;
1426                         sys->molblock[mb + 1].nmol -= 1;
1427                     }
1428
1429                     /* Create a copy of a moltype for a molecule
1430                      * containing QM atoms and append it in the end of the list
1431                      */
1432                     std::vector<gmx_moltype_t> temp(sys->moltype.size());
1433                     for (size_t i = 0; i < sys->moltype.size(); ++i)
1434                     {
1435                         copy_moltype(&sys->moltype[i], &temp[i]);
1436                     }
1437                     sys->moltype.resize(sys->moltype.size() + 1);
1438                     for (size_t i = 0; i < temp.size(); ++i)
1439                     {
1440                         copy_moltype(&temp[i], &sys->moltype[i]);
1441                     }
1442                     copy_moltype(&sys->moltype[molb->type], &sys->moltype.back());
1443                     /* Copy the exclusions to a new array, since this is the only
1444                      * thing that needs to be modified for QMMM.
1445                      */
1446                     sys->moltype.back().excls = sys->moltype[molb->type].excls;
1447                     /* Set the molecule type for the QMMM molblock */
1448                     molb->type = sys->moltype.size() - 1;
1449                 }
1450                 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, logger);
1451             }
1452             if (grpnr)
1453             {
1454                 grpnr += nat_mol;
1455             }
1456             index_offset += nat_mol;
1457         }
1458     }
1459 }