src/gromacs/gmxpreprocess/topio.c

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