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