src/gromacs/gmxpreprocess/topio.cpp

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