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