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