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