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