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