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