src/gromacs/gmxpreprocess/topio.c

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