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