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