src/gromacs/gmxpreprocess/topio.cpp

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