src/gromacs/gmxpreprocess/topio.cpp

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