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