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