src/gromacs/gmxpreprocess/topio.cpp

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