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