src/gromacs/gmxpreprocess/topio.cpp

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