src/gromacs/gmxpreprocess/toppush.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 <ctype.h>
  41 #include <math.h>
  42
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "string2.h"
  47 #include "names.h"
  48 #include "toputil.h"
  49 #include "toppush.h"
  50 #include "topdirs.h"
  51 #include "readir.h"
  52 #include "symtab.h"
  53 #include "gmx_fatal.h"
  54 #include "warninp.h"
  55 #include "gpp_atomtype.h"
  56 #include "gpp_bond_atomtype.h"
  57
  58 void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype,
  59                        warninp_t wi)
  60 {
  61     int   i, j, k = -1, nf;
  62     int   nr, nrfp;
  63     real  c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2;
  64     char  errbuf[256];
  65
  66     /* Lean mean shortcuts */
  67     nr   = get_atomtype_ntypes(atype);
  68     nrfp = NRFP(ftype);
  69     snew(plist->param, nr*nr);
  70     plist->nr = nr*nr;
  71
  72     /* Fill the matrix with force parameters */
  73     switch (ftype)
  74     {
  75         case F_LJ:
  76             switch (comb)
  77             {
  78                 case eCOMB_GEOMETRIC:
  79                     /* Gromos rules */
  80                     for (i = k = 0; (i < nr); i++)
  81                     {
  82                         for (j = 0; (j < nr); j++, k++)
  83                         {
  84                             for (nf = 0; (nf < nrfp); nf++)
  85                             {
  86                                 ci = get_atomtype_nbparam(i, nf, atype);
  87                                 cj = get_atomtype_nbparam(j, nf, atype);
  88                                 c  = sqrt(ci * cj);
  89                                 plist->param[k].c[nf]      = c;
  90                             }
  91                         }
  92                     }
  93                     break;
  94
  95                 case eCOMB_ARITHMETIC:
  96                     /* c0 and c1 are epsilon and sigma */
  97                     for (i = k = 0; (i < nr); i++)
  98                     {
  99                         for (j = 0; (j < nr); j++, k++)
 100                         {
 101                             ci0                  = get_atomtype_nbparam(i, 0, atype);
 102                             cj0                  = get_atomtype_nbparam(j, 0, atype);
 103                             ci1                  = get_atomtype_nbparam(i, 1, atype);
 104                             cj1                  = get_atomtype_nbparam(j, 1, atype);
 105                             plist->param[k].c[0] = (ci0+cj0)*0.5;
 106                             plist->param[k].c[1] = sqrt(ci1*cj1);
 107                         }
 108                     }
 109
 110                     break;
 111                 case eCOMB_GEOM_SIG_EPS:
 112                     /* c0 and c1 are epsilon and sigma */
 113                     for (i = k = 0; (i < nr); i++)
 114                     {
 115                         for (j = 0; (j < nr); j++, k++)
 116                         {
 117                             ci0                  = get_atomtype_nbparam(i, 0, atype);
 118                             cj0                  = get_atomtype_nbparam(j, 0, atype);
 119                             ci1                  = get_atomtype_nbparam(i, 1, atype);
 120                             cj1                  = get_atomtype_nbparam(j, 1, atype);
 121                             plist->param[k].c[0] = sqrt(ci0*cj0);
 122                             plist->param[k].c[1] = sqrt(ci1*cj1);
 123                         }
 124                     }
 125
 126                     break;
 127                 default:
 128                     gmx_fatal(FARGS, "No such combination rule %d", comb);
 129             }
 130             if (plist->nr != k)
 131             {
 132                 gmx_incons("Topology processing, generate nb parameters");
 133             }
 134             break;
 135
 136         case F_BHAM:
 137             /* Buckingham rules */
 138             for (i = k = 0; (i < nr); i++)
 139             {
 140                 for (j = 0; (j < nr); j++, k++)
 141                 {
 142                     ci0                  = get_atomtype_nbparam(i, 0, atype);
 143                     cj0                  = get_atomtype_nbparam(j, 0, atype);
 144                     ci2                  = get_atomtype_nbparam(i, 2, atype);
 145                     cj2                  = get_atomtype_nbparam(j, 2, atype);
 146                     bi                   = get_atomtype_nbparam(i, 1, atype);
 147                     bj                   = get_atomtype_nbparam(j, 1, atype);
 148                     plist->param[k].c[0] = sqrt(ci0 * cj0);
 149                     if ((bi == 0) || (bj == 0))
 150                     {
 151                         plist->param[k].c[1] = 0;
 152                     }
 153                     else
 154                     {
 155                         plist->param[k].c[1] = 2.0/(1/bi+1/bj);
 156                     }
 157                     plist->param[k].c[2] = sqrt(ci2 * cj2);
 158                 }
 159             }
 160
 161             break;
 162         default:
 163             sprintf(errbuf, "Invalid nonbonded type %s",
 164                     interaction_function[ftype].longname);
 165             warning_error(wi, errbuf);
 166     }
 167 }
 168
 169 static void realloc_nb_params(gpp_atomtype_t at,
 170                               t_nbparam ***nbparam, t_nbparam ***pair)
 171 {
 172     /* Add space in the non-bonded parameters matrix */
 173     int atnr = get_atomtype_ntypes(at);
 174     srenew(*nbparam, atnr);
 175     snew((*nbparam)[atnr-1], atnr);
 176     if (pair)
 177     {
 178         srenew(*pair, atnr);
 179         snew((*pair)[atnr-1], atnr);
 180     }
 181 }
 182
 183 static void copy_B_from_A(int ftype, double *c)
 184 {
 185     int nrfpA, nrfpB, i;
 186
 187     nrfpA = NRFPA(ftype);
 188     nrfpB = NRFPB(ftype);
 189
 190     /* Copy the B parameters from the first nrfpB A parameters */
 191     for (i = 0; (i < nrfpB); i++)
 192     {
 193         c[nrfpA+i] = c[i];
 194     }
 195 }
 196
 197 void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat,
 198               char *line, int nb_funct,
 199               t_nbparam ***nbparam, t_nbparam ***pair,
 200               warninp_t wi)
 201 {
 202     typedef struct {
 203         const char *entry;
 204         int         ptype;
 205     } t_xlate;
 206     t_xlate    xl[eptNR] = {
 207         { "A",   eptAtom },
 208         { "N",   eptNucleus },
 209         { "S",   eptShell },
 210         { "B",   eptBond },
 211         { "V",   eptVSite },
 212     };
 213
 214     int        nr, i, nfields, j, pt, nfp0 = -1;
 215     int        batype_nr, nread;
 216     char       type[STRLEN], btype[STRLEN], ptype[STRLEN];
 217     double     m, q;
 218     double     c[MAXFORCEPARAM];
 219     double     radius, vol, surftens, gb_radius, S_hct;
 220     char       tmpfield[12][100]; /* Max 12 fields of width 100 */
 221     char       errbuf[256];
 222     t_atom    *atom;
 223     t_param   *param;
 224     int        atomnr;
 225     gmx_bool   have_atomic_number;
 226     gmx_bool   have_bonded_type;
 227
 228     snew(atom, 1);
 229     snew(param, 1);
 230
 231     /* First assign input line to temporary array */
 232     nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s",
 233                      tmpfield[0], tmpfield[1], tmpfield[2], tmpfield[3], tmpfield[4], tmpfield[5],
 234                      tmpfield[6], tmpfield[7], tmpfield[8], tmpfield[9], tmpfield[10], tmpfield[11]);
 235
 236     /* Comments on optional fields in the atomtypes section:
 237      *
 238      * The force field format is getting a bit old. For OPLS-AA we needed
 239      * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
 240      * we also needed the atomic numbers.
 241      * To avoid making all old or user-generated force fields unusable we
 242      * have introduced both these quantities as optional columns, and do some
 243      * acrobatics to check whether they are present or not.
 244      * This will all look much nicer when we switch to XML... sigh.
 245      *
 246      * Field 0 (mandatory) is the nonbonded type name. (string)
 247      * Field 1 (optional)  is the bonded type (string)
 248      * Field 2 (optional)  is the atomic number (int)
 249      * Field 3 (mandatory) is the mass (numerical)
 250      * Field 4 (mandatory) is the charge (numerical)
 251      * Field 5 (mandatory) is the particle type (single character)
 252      * This is followed by a number of nonbonded parameters.
 253      *
 254      * The safest way to identify the format is the particle type field.
 255      *
 256      * So, here is what we do:
 257      *
 258      * A. Read in the first six fields as strings
 259      * B. If field 3 (starting from 0) is a single char, we have neither
 260      *    bonded_type or atomic numbers.
 261      * C. If field 5 is a single char we have both.
 262      * D. If field 4 is a single char we check field 1. If this begins with
 263      *    an alphabetical character we have bonded types, otherwise atomic numbers.
 264      */
 265
 266     if (nfields < 6)
 267     {
 268         too_few(wi);
 269         return;
 270     }
 271
 272     if ( (strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0]) )
 273     {
 274         have_bonded_type   = TRUE;
 275         have_atomic_number = TRUE;
 276     }
 277     else if ( (strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0]) )
 278     {
 279         have_bonded_type   = FALSE;
 280         have_atomic_number = FALSE;
 281     }
 282     else
 283     {
 284         have_bonded_type   = ( isalpha(tmpfield[1][0]) != 0 );
 285         have_atomic_number = !have_bonded_type;
 286     }
 287
 288     /* optional fields */
 289     surftens  = -1;
 290     vol       = -1;
 291     radius    = -1;
 292     gb_radius = -1;
 293     atomnr    = -1;
 294     S_hct     = -1;
 295
 296     switch (nb_funct)
 297     {
 298
 299         case F_LJ:
 300             nfp0 = 2;
 301
 302             if (have_atomic_number)
 303             {
 304                 if (have_bonded_type)
 305                 {
 306                     nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
 307                                    type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1],
 308                                    &radius, &vol, &surftens, &gb_radius);
 309                     if (nread < 8)
 310                     {
 311                         too_few(wi);
 312                         return;
 313                     }
 314                 }
 315                 else
 316                 {
 317                     /* have_atomic_number && !have_bonded_type */
 318                     nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
 319                                    type, &atomnr, &m, &q, ptype, &c[0], &c[1],
 320                                    &radius, &vol, &surftens, &gb_radius);
 321                     if (nread < 7)
 322                     {
 323                         too_few(wi);
 324                         return;
 325                     }
 326                 }
 327             }
 328             else
 329             {
 330                 if (have_bonded_type)
 331                 {
 332                     /* !have_atomic_number && have_bonded_type */
 333                     nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
 334                                    type, btype, &m, &q, ptype, &c[0], &c[1],
 335                                    &radius, &vol, &surftens, &gb_radius);
 336                     if (nread < 7)
 337                     {
 338                         too_few(wi);
 339                         return;
 340                     }
 341                 }
 342                 else
 343                 {
 344                     /* !have_atomic_number && !have_bonded_type */
 345                     nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
 346                                    type, &m, &q, ptype, &c[0], &c[1],
 347                                    &radius, &vol, &surftens, &gb_radius);
 348                     if (nread < 6)
 349                     {
 350                         too_few(wi);
 351                         return;
 352                     }
 353                 }
 354             }
 355
 356             if (!have_bonded_type)
 357             {
 358                 strcpy(btype, type);
 359             }
 360
 361             if (!have_atomic_number)
 362             {
 363                 atomnr = -1;
 364             }
 365
 366             break;
 367
 368         case F_BHAM:
 369             nfp0 = 3;
 370
 371             if (have_atomic_number)
 372             {
 373                 if (have_bonded_type)
 374                 {
 375                     nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
 376                                    type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2],
 377                                    &radius, &vol, &surftens, &gb_radius);
 378                     if (nread < 9)
 379                     {
 380                         too_few(wi);
 381                         return;
 382                     }
 383                 }
 384                 else
 385                 {
 386                     /* have_atomic_number && !have_bonded_type */
 387                     nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
 388                                    type, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2],
 389                                    &radius, &vol, &surftens, &gb_radius);
 390                     if (nread < 8)
 391                     {
 392                         too_few(wi);
 393                         return;
 394                     }
 395                 }
 396             }
 397             else
 398             {
 399                 if (have_bonded_type)
 400                 {
 401                     /* !have_atomic_number && have_bonded_type */
 402                     nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
 403                                    type, btype, &m, &q, ptype, &c[0], &c[1], &c[2],
 404                                    &radius, &vol, &surftens, &gb_radius);
 405                     if (nread < 8)
 406                     {
 407                         too_few(wi);
 408                         return;
 409                     }
 410                 }
 411                 else
 412                 {
 413                     /* !have_atomic_number && !have_bonded_type */
 414                     nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
 415                                    type, &m, &q, ptype, &c[0], &c[1], &c[2],
 416                                    &radius, &vol, &surftens, &gb_radius);
 417                     if (nread < 7)
 418                     {
 419                         too_few(wi);
 420                         return;
 421                     }
 422                 }
 423             }
 424
 425             if (!have_bonded_type)
 426             {
 427                 strcpy(btype, type);
 428             }
 429
 430             if (!have_atomic_number)
 431             {
 432                 atomnr = -1;
 433             }
 434
 435             break;
 436
 437         default:
 438             gmx_fatal(FARGS, "Invalid function type %d in push_at %s %d", nb_funct,
 439                       __FILE__, __LINE__);
 440     }
 441     for (j = nfp0; (j < MAXFORCEPARAM); j++)
 442     {
 443         c[j] = 0.0;
 444     }
 445
 446     if (strlen(type) == 1 && isdigit(type[0]))
 447     {
 448         gmx_fatal(FARGS, "Atom type names can't be single digits.");
 449     }
 450
 451     if (strlen(btype) == 1 && isdigit(btype[0]))
 452     {
 453         gmx_fatal(FARGS, "Bond atom type names can't be single digits.");
 454     }
 455
 456     /* Hack to read old topologies */
 457     if (gmx_strcasecmp(ptype, "D") == 0)
 458     {
 459         sprintf(ptype, "V");
 460     }
 461     for (j = 0; (j < eptNR); j++)
 462     {
 463         if (gmx_strcasecmp(ptype, xl[j].entry) == 0)
 464         {
 465             break;
 466         }
 467     }
 468     if (j == eptNR)
 469     {
 470         gmx_fatal(FARGS, "Invalid particle type %s on line %s",
 471                   ptype, line);
 472     }
 473     pt = xl[j].ptype;
 474     if (debug)
 475     {
 476         fprintf(debug, "ptype: %s\n", ptype_str[pt]);
 477     }
 478
 479     atom->q     = q;
 480     atom->m     = m;
 481     atom->ptype = pt;
 482     for (i = 0; (i < MAXFORCEPARAM); i++)
 483     {
 484         param->c[i] = c[i];
 485     }
 486
 487     if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET)
 488     {
 489         add_bond_atomtype(bat, symtab, btype);
 490     }
 491     batype_nr = get_bond_atomtype_type(btype, bat);
 492
 493     if ((nr = get_atomtype_type(type, at)) != NOTSET)
 494     {
 495         sprintf(errbuf, "Overriding atomtype %s", type);
 496         warning(wi, errbuf);
 497         if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr,
 498                                radius, vol, surftens, atomnr, gb_radius, S_hct)) == NOTSET)
 499         {
 500             gmx_fatal(FARGS, "Replacing atomtype %s failed", type);
 501         }
 502     }
 503     else if ((nr = add_atomtype(at, symtab, atom, type, param,
 504                                 batype_nr, radius, vol,
 505                                 surftens, atomnr, gb_radius, S_hct)) == NOTSET)
 506     {
 507         gmx_fatal(FARGS, "Adding atomtype %s failed", type);
 508     }
 509     else
 510     {
 511         /* Add space in the non-bonded parameters matrix */
 512         realloc_nb_params(at, nbparam, pair);
 513     }
 514     sfree(atom);
 515     sfree(param);
 516 }
 517
 518 static void push_bondtype(t_params     *       bt,
 519                           t_param     *        b,
 520                           int                  nral,
 521                           int                  ftype,
 522                           gmx_bool             bAllowRepeat,
 523                           char     *           line,
 524                           warninp_t            wi)
 525 {
 526     int      i, j;
 527     gmx_bool bTest, bFound, bCont, bId;
 528     int      nr   = bt->nr;
 529     int      nrfp = NRFP(ftype);
 530     char     errbuf[256];
 531
 532     /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
 533        are on directly _adjacent_ lines.
 534      */
 535
 536     /* First check if our atomtypes are _identical_ (not reversed) to the previous
 537        entry. If they are not identical we search for earlier duplicates. If they are
 538        we can skip it, since we already searched for the first line
 539        in this group.
 540      */
 541
 542     bFound = FALSE;
 543     bCont  = FALSE;
 544
 545     if (bAllowRepeat && nr > 1)
 546     {
 547         for (j = 0, bCont = TRUE; (j < nral); j++)
 548         {
 549             bCont = bCont && (b->a[j] == bt->param[nr-2].a[j]);
 550         }
 551     }
 552
 553     /* Search for earlier duplicates if this entry was not a continuation
 554        from the previous line.
 555      */
 556     if (!bCont)
 557     {
 558         bFound = FALSE;
 559         for (i = 0; (i < nr); i++)
 560         {
 561             bTest = TRUE;
 562             for (j = 0; (j < nral); j++)
 563             {
 564                 bTest = (bTest && (b->a[j] == bt->param[i].a[j]));
 565             }
 566             if (!bTest)
 567             {
 568                 bTest = TRUE;
 569                 for (j = 0; (j < nral); j++)
 570                 {
 571                     bTest = (bTest && (b->a[nral-1-j] == bt->param[i].a[j]));
 572                 }
 573             }
 574             if (bTest)
 575             {
 576                 if (!bFound)
 577                 {
 578                     bId = TRUE;
 579                     for (j = 0; (j < nrfp); j++)
 580                     {
 581                         bId = bId && (bt->param[i].c[j] == b->c[j]);
 582                     }
 583                     if (!bId)
 584                     {
 585                         sprintf(errbuf, "Overriding %s parameters.%s",
 586                                 interaction_function[ftype].longname,
 587                                 (ftype == F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : "");
 588                         warning(wi, errbuf);
 589                         fprintf(stderr, "  old:");
 590                         for (j = 0; (j < nrfp); j++)
 591                         {
 592                             fprintf(stderr, " %g", bt->param[i].c[j]);
 593                         }
 594                         fprintf(stderr, " \n  new: %s\n\n", line);
 595                     }
 596                 }
 597                 /* Overwrite it! */
 598                 for (j = 0; (j < nrfp); j++)
 599                 {
 600                     bt->param[i].c[j] = b->c[j];
 601                 }
 602                 bFound = TRUE;
 603             }
 604         }
 605     }
 606     if (!bFound)
 607     {
 608         /* alloc */
 609         pr_alloc (2, bt);
 610
 611         /* fill the arrays up and down */
 612         memcpy(bt->param[bt->nr].c,  b->c, sizeof(b->c));
 613         memcpy(bt->param[bt->nr].a,  b->a, sizeof(b->a));
 614         memcpy(bt->param[bt->nr+1].c, b->c, sizeof(b->c));
 615
 616         /* The definitions of linear angles depend on the order of atoms,
 617          * that means that for atoms i-j-k, with certain parameter a, the
 618          * corresponding k-j-i angle will have parameter 1-a.
 619          */
 620         if (ftype == F_LINEAR_ANGLES)
 621         {
 622             bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0];
 623             bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2];
 624         }
 625
 626         for (j = 0; (j < nral); j++)
 627         {
 628             bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
 629         }
 630
 631         bt->nr += 2;
 632     }
 633 }
 634
 635 void push_bt(directive d, t_params bt[], int nral,
 636              gpp_atomtype_t at,
 637              t_bond_atomtype bat, char *line,
 638              warninp_t wi)
 639 {
 640     const char *formal[MAXATOMLIST+1] = {
 641         "%s",
 642         "%s%s",
 643         "%s%s%s",
 644         "%s%s%s%s",
 645         "%s%s%s%s%s",
 646         "%s%s%s%s%s%s",
 647         "%s%s%s%s%s%s%s"
 648     };
 649     const char *formnl[MAXATOMLIST+1] = {
 650         "%*s",
 651         "%*s%*s",
 652         "%*s%*s%*s",
 653         "%*s%*s%*s%*s",
 654         "%*s%*s%*s%*s%*s",
 655         "%*s%*s%*s%*s%*s%*s",
 656         "%*s%*s%*s%*s%*s%*s%*s"
 657     };
 658     const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
 659     int         i, ft, ftype, nn, nrfp, nrfpA, nrfpB;
 660     char        f1[STRLEN];
 661     char        alc[MAXATOMLIST+1][20];
 662     /* One force parameter more, so we can check if we read too many */
 663     double      c[MAXFORCEPARAM+1];
 664     t_param     p;
 665     char        errbuf[256];
 666
 667     if ((bat && at) || (!bat && !at))
 668     {
 669         gmx_incons("You should pass either bat or at to push_bt");
 670     }
 671
 672     /* Make format string (nral ints+functype) */
 673     if ((nn = sscanf(line, formal[nral],
 674                      alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral+1)
 675     {
 676         sprintf(errbuf, "Not enough atomtypes (%d instead of %d)", nn-1, nral);
 677         warning_error(wi, errbuf);
 678         return;
 679     }
 680
 681     ft    = strtol(alc[nral], NULL, 10);
 682     ftype = ifunc_index(d, ft);
 683     nrfp  = NRFP(ftype);
 684     nrfpA = interaction_function[ftype].nrfpA;
 685     nrfpB = interaction_function[ftype].nrfpB;
 686     strcpy(f1, formnl[nral]);
 687     strcat(f1, formlf);
 688     if ((nn = sscanf(line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9], &c[10], &c[11], &c[12]))
 689         != nrfp)
 690     {
 691         if (nn == nrfpA)
 692         {
 693             /* Copy the B-state from the A-state */
 694             copy_B_from_A(ftype, c);
 695         }
 696         else
 697         {
 698             if (nn < nrfpA)
 699             {
 700                 warning_error(wi, "Not enough parameters");
 701             }
 702             else if (nn > nrfpA && nn < nrfp)
 703             {
 704                 warning_error(wi, "Too many parameters or not enough parameters for topology B");
 705             }
 706             else if (nn > nrfp)
 707             {
 708                 warning_error(wi, "Too many parameters");
 709             }
 710             for (i = nn; (i < nrfp); i++)
 711             {
 712                 c[i] = 0.0;
 713             }
 714         }
 715     }
 716     for (i = 0; (i < nral); i++)
 717     {
 718         if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET))
 719         {
 720             gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]);
 721         }
 722         else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
 723         {
 724             gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]);
 725         }
 726     }
 727     for (i = 0; (i < nrfp); i++)
 728     {
 729         p.c[i] = c[i];
 730     }
 731     push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
 732 }
 733
 734
 735 void push_dihedraltype(directive d, t_params bt[],
 736                        t_bond_atomtype bat, char *line,
 737                        warninp_t wi)
 738 {
 739     const char  *formal[MAXATOMLIST+1] = {
 740         "%s",
 741         "%s%s",
 742         "%s%s%s",
 743         "%s%s%s%s",
 744         "%s%s%s%s%s",
 745         "%s%s%s%s%s%s",
 746         "%s%s%s%s%s%s%s"
 747     };
 748     const char  *formnl[MAXATOMLIST+1] = {
 749         "%*s",
 750         "%*s%*s",
 751         "%*s%*s%*s",
 752         "%*s%*s%*s%*s",
 753         "%*s%*s%*s%*s%*s",
 754         "%*s%*s%*s%*s%*s%*s",
 755         "%*s%*s%*s%*s%*s%*s%*s"
 756     };
 757     const char  *formlf[MAXFORCEPARAM] = {
 758         "%lf",
 759         "%lf%lf",
 760         "%lf%lf%lf",
 761         "%lf%lf%lf%lf",
 762         "%lf%lf%lf%lf%lf",
 763         "%lf%lf%lf%lf%lf%lf",
 764         "%lf%lf%lf%lf%lf%lf%lf",
 765         "%lf%lf%lf%lf%lf%lf%lf%lf",
 766         "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 767         "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 768         "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 769         "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
 770     };
 771     int          i, ft, ftype, nn, nrfp, nrfpA, nrfpB, nral;
 772     char         f1[STRLEN];
 773     char         alc[MAXATOMLIST+1][20];
 774     double       c[MAXFORCEPARAM];
 775     t_param      p;
 776     gmx_bool     bAllowRepeat;
 777     char         errbuf[256];
 778
 779     /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
 780      *
 781      * We first check for 2 atoms with the 3th column being an integer
 782      * defining the type. If this isn't the case, we try it with 4 atoms
 783      * and the 5th column defining the dihedral type.
 784      */
 785     nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]);
 786     if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0]))
 787     {
 788         nral  = 2;
 789         ft    = strtol(alc[nral], NULL, 10);
 790         /* Move atom types around a bit and use 'X' for wildcard atoms
 791          * to create a 4-atom dihedral definition with arbitrary atoms in
 792          * position 1 and 4.
 793          */
 794         if (alc[2][0] == '2')
 795         {
 796             /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
 797             strcpy(alc[3], alc[1]);
 798             sprintf(alc[2], "X");
 799             sprintf(alc[1], "X");
 800             /* alc[0] stays put */
 801         }
 802         else
 803         {
 804             /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
 805             sprintf(alc[3], "X");
 806             strcpy(alc[2], alc[1]);
 807             strcpy(alc[1], alc[0]);
 808             sprintf(alc[0], "X");
 809         }
 810     }
 811     else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0]))
 812     {
 813         nral  = 4;
 814         ft    = strtol(alc[nral], NULL, 10);
 815     }
 816     else
 817     {
 818         sprintf(errbuf, "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn-1);
 819         warning_error(wi, errbuf);
 820         return;
 821     }
 822
 823     if (ft == 9)
 824     {
 825         /* Previously, we have always overwritten parameters if e.g. a torsion
 826            with the same atomtypes occurs on multiple lines. However, CHARMM and
 827            some other force fields specify multiple dihedrals over some bonds,
 828            including cosines with multiplicity 6 and somethimes even higher.
 829            Thus, they cannot be represented with Ryckaert-Bellemans terms.
 830            To add support for these force fields, Dihedral type 9 is identical to
 831            normal proper dihedrals, but repeated entries are allowed.
 832          */
 833         bAllowRepeat = TRUE;
 834         ft           = 1;
 835     }
 836     else
 837     {
 838         bAllowRepeat = FALSE;
 839     }
 840
 841
 842     ftype = ifunc_index(d, ft);
 843     nrfp  = NRFP(ftype);
 844     nrfpA = interaction_function[ftype].nrfpA;
 845     nrfpB = interaction_function[ftype].nrfpB;
 846
 847     strcpy(f1, formnl[nral]);
 848     strcat(f1, formlf[nrfp-1]);
 849
 850     /* Check number of parameters given */
 851     if ((nn = sscanf(line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9], &c[10], &c[11]))
 852         != nrfp)
 853     {
 854         if (nn == nrfpA)
 855         {
 856             /* Copy the B-state from the A-state */
 857             copy_B_from_A(ftype, c);
 858         }
 859         else
 860         {
 861             if (nn < nrfpA)
 862             {
 863                 warning_error(wi, "Not enough parameters");
 864             }
 865             else if (nn > nrfpA && nn < nrfp)
 866             {
 867                 warning_error(wi, "Too many parameters or not enough parameters for topology B");
 868             }
 869             else if (nn > nrfp)
 870             {
 871                 warning_error(wi, "Too many parameters");
 872             }
 873             for (i = nn; (i < nrfp); i++)
 874             {
 875                 c[i] = 0.0;
 876             }
 877         }
 878     }
 879
 880     for (i = 0; (i < 4); i++)
 881     {
 882         if (!strcmp(alc[i], "X"))
 883         {
 884             p.a[i] = -1;
 885         }
 886         else
 887         {
 888             if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)
 889             {
 890                 gmx_fatal(FARGS, "Unknown bond_atomtype %s", alc[i]);
 891             }
 892         }
 893     }
 894     for (i = 0; (i < nrfp); i++)
 895     {
 896         p.c[i] = c[i];
 897     }
 898     /* Always use 4 atoms here, since we created two wildcard atoms
 899      * if there wasn't of them 4 already.
 900      */
 901     push_bondtype (&(bt[ftype]), &p, 4, ftype, bAllowRepeat, line, wi);
 902 }
 903
 904
 905 void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype,
 906               char *pline, int nb_funct,
 907               warninp_t wi)
 908 {
 909     /* swap the atoms */
 910     const char *form2 = "%*s%*s%*s%lf%lf";
 911     const char *form3 = "%*s%*s%*s%lf%lf%lf";
 912     const char *form4 = "%*s%*s%*s%lf%lf%lf%lf";
 913     const char *form5 = "%*s%*s%*s%lf%lf%lf%lf%lf";
 914     char        a0[80], a1[80];
 915     int         i, f, n, ftype, atnr, nrfp;
 916     double      c[4], dum;
 917     real        cr[4], sig6;
 918     atom_id     ai, aj;
 919     t_nbparam  *nbp;
 920     gmx_bool    bId;
 921     char        errbuf[256];
 922
 923     if (sscanf (pline, "%s%s%d", a0, a1, &f) != 3)
 924     {
 925         too_few(wi);
 926         return;
 927     }
 928
 929     ftype = ifunc_index(d, f);
 930
 931     if (ftype != nb_funct)
 932     {
 933         sprintf(errbuf, "Trying to add %s while the default nonbond type is %s",
 934                 interaction_function[ftype].longname,
 935                 interaction_function[nb_funct].longname);
 936         warning_error(wi, errbuf);
 937         return;
 938     }
 939
 940     /* Get the force parameters */
 941     nrfp = NRFP(ftype);
 942     if (ftype == F_LJ14)
 943     {
 944         n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]);
 945         if (n < 2)
 946         {
 947             too_few(wi);
 948             return;
 949         }
 950         /* When the B topology parameters are not set,
 951          * copy them from topology A
 952          */
 953         for (i = n; i < nrfp; i++)
 954         {
 955             c[i] = c[i-2];
 956         }
 957     }
 958     else if (ftype == F_LJC14_Q)
 959     {
 960         n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum);
 961         if (n != 4)
 962         {
 963             incorrect_n_param(wi);
 964             return;
 965         }
 966     }
 967     else if (nrfp == 2)
 968     {
 969         if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2)
 970         {
 971             incorrect_n_param(wi);
 972             return;
 973         }
 974     }
 975     else if (nrfp == 3)
 976     {
 977         if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3)
 978         {
 979             incorrect_n_param(wi);
 980             return;
 981         }
 982     }
 983     else
 984     {
 985         gmx_fatal(FARGS, "Number of force parameters for nonbonded interactions is %d"
 986                   " in file %s, line %d", nrfp, __FILE__, __LINE__);
 987     }
 988     for (i = 0; (i < nrfp); i++)
 989     {
 990         cr[i] = c[i];
 991     }
 992
 993     /* Put the parameters in the matrix */
 994     if ((ai = get_atomtype_type (a0, atype)) == NOTSET)
 995     {
 996         gmx_fatal(FARGS, "Atomtype %s not found", a0);
 997     }
 998     if ((aj = get_atomtype_type (a1, atype)) == NOTSET)
 999     {
1000         gmx_fatal(FARGS, "Atomtype %s not found", a1);
1001     }
1002     nbp = &(nbt[max(ai, aj)][min(ai, aj)]);
1003
1004     if (nbp->bSet)
1005     {
1006         bId = TRUE;
1007         for (i = 0; i < nrfp; i++)
1008         {
1009             bId = bId && (nbp->c[i] == cr[i]);
1010         }
1011         if (!bId)
1012         {
1013             sprintf(errbuf, "Overriding non-bonded parameters,");
1014             warning(wi, errbuf);
1015             fprintf(stderr, "  old:");
1016             for (i = 0; i < nrfp; i++)
1017             {
1018                 fprintf(stderr, " %g", nbp->c[i]);
1019             }
1020             fprintf(stderr, " new\n%s\n", pline);
1021         }
1022     }
1023     nbp->bSet = TRUE;
1024     for (i = 0; i < nrfp; i++)
1025     {
1026         nbp->c[i] = cr[i];
1027     }
1028 }
1029
1030 void
1031 push_gb_params (gpp_atomtype_t at, char *line,
1032                 warninp_t wi)
1033 {
1034     int    nfield;
1035     int    atype;
1036     double radius, vol, surftens, gb_radius, S_hct;
1037     char   atypename[STRLEN];
1038     char   errbuf[STRLEN];
1039
1040     if ( (nfield = sscanf(line, "%s%lf%lf%lf%lf%lf", atypename, &radius, &vol, &surftens, &gb_radius, &S_hct)) != 6)
1041     {
1042         sprintf(errbuf, "Too few gb parameters for type %s\n", atypename);
1043         warning(wi, errbuf);
1044     }
1045
1046     /* Search for atomtype */
1047     atype = get_atomtype_type(atypename, at);
1048
1049     if (atype == NOTSET)
1050     {
1051         printf("Couldn't find topology match for atomtype %s\n", atypename);
1052         abort();
1053     }
1054
1055     set_atomtype_gbparam(at, atype, radius, vol, surftens, gb_radius, S_hct);
1056 }
1057
1058 void
1059 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
1060               t_bond_atomtype bat, char *line,
1061               warninp_t wi)
1062 {
1063     const char  *formal = "%s%s%s%s%s%s%s%s";
1064
1065     int          i, j, ft, ftype, nn, nrfp, nrfpA, nrfpB;
1066     int          start;
1067     int          nxcmap, nycmap, ncmap, read_cmap, sl, nct;
1068     char         s[20], alc[MAXATOMLIST+2][20];
1069     t_param      p;
1070     gmx_bool     bAllowRepeat;
1071     char         errbuf[256];
1072
1073     /* Keep the compiler happy */
1074     read_cmap = 0;
1075     start     = 0;
1076
1077     if ((nn = sscanf(line, formal, alc[0], alc[1], alc[2], alc[3], alc[4], alc[5], alc[6], alc[7])) != nral+3)
1078     {
1079         sprintf(errbuf, "Incorrect number of atomtypes for cmap (%d instead of 5)", nn-1);
1080         warning_error(wi, errbuf);
1081         return;
1082     }
1083
1084     /* Compute an offset for each line where the cmap parameters start
1085      * ie. where the atom types and grid spacing information ends
1086      */
1087     for (i = 0; i < nn; i++)
1088     {
1089         start += (int)strlen(alc[i]);
1090     }
1091
1092     /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */
1093     /* start is the position on the line where we start to read the actual cmap grid data from the itp file */
1094     start = start + nn -1;
1095
1096     ft     = strtol(alc[nral], NULL, 10);
1097     nxcmap = strtol(alc[nral+1], NULL, 10);
1098     nycmap = strtol(alc[nral+2], NULL, 10);
1099
1100     /* Check for equal grid spacing in x and y dims */
1101     if (nxcmap != nycmap)
1102     {
1103         gmx_fatal(FARGS, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
1104     }
1105
1106     ncmap  = nxcmap*nycmap;
1107     ftype  = ifunc_index(d, ft);
1108     nrfpA  = strtol(alc[6], NULL, 10)*strtol(alc[6], NULL, 10);
1109     nrfpB  = strtol(alc[7], NULL, 10)*strtol(alc[7], NULL, 10);
1110     nrfp   = nrfpA+nrfpB;
1111
1112     /* Allocate memory for the CMAP grid */
1113     bt->ncmap += nrfp;
1114     srenew(bt->cmap, bt->ncmap);
1115
1116     /* Read in CMAP parameters */
1117     sl = 0;
1118     for (i = 0; i < ncmap; i++)
1119     {
1120         while (isspace(*(line+start+sl)))
1121         {
1122             sl++;
1123         }
1124         nn  = sscanf(line+start+sl, " %s ", s);
1125         sl += strlen(s);
1126         bt->cmap[i+(bt->ncmap)-nrfp] = strtod(s, NULL);
1127
1128         if (nn == 1)
1129         {
1130             read_cmap++;
1131         }
1132         else
1133         {
1134             gmx_fatal(FARGS, "Error in reading cmap parameter for angle %s %s %s %s %s", alc[0], alc[1], alc[2], alc[3], alc[4]);
1135         }
1136
1137     }
1138
1139     /* Check do that we got the number of parameters we expected */
1140     if (read_cmap == nrfpA)
1141     {
1142         for (i = 0; i < ncmap; i++)
1143         {
1144             bt->cmap[i+ncmap] = bt->cmap[i];
1145         }
1146     }
1147     else
1148     {
1149         if (read_cmap < nrfpA)
1150         {
1151             warning_error(wi, "Not enough cmap parameters");
1152         }
1153         else if (read_cmap > nrfpA && read_cmap < nrfp)
1154         {
1155             warning_error(wi, "Too many cmap parameters or not enough parameters for topology B");
1156         }
1157         else if (read_cmap > nrfp)
1158         {
1159             warning_error(wi, "Too many cmap parameters");
1160         }
1161     }
1162
1163
1164     /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1165      * so we can safely assign them each time
1166      */
1167     bt->grid_spacing = nxcmap;     /* Or nycmap, they need to be equal */
1168     bt->nc           = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1169     nct              = (nral+1) * bt->nc;
1170
1171     /* Allocate memory for the cmap_types information */
1172     srenew(bt->cmap_types, nct);
1173
1174     for (i = 0; (i < nral); i++)
1175     {
1176         if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1177         {
1178             gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]);
1179         }
1180         else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1181         {
1182             gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]);
1183         }
1184
1185         /* Assign a grid number to each cmap_type */
1186         bt->cmap_types[bt->nct++] = get_bond_atomtype_type(alc[i], bat);
1187     }
1188
1189     /* Assign a type number to this cmap */
1190     bt->cmap_types[bt->nct++] = bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1191
1192     /* Check for the correct number of atoms (again) */
1193     if (bt->nct != nct)
1194     {
1195         gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt->nc);
1196     }
1197
1198     /* Is this correct?? */
1199     for (i = 0; i < MAXFORCEPARAM; i++)
1200     {
1201         p.c[i] = NOTSET;
1202     }
1203
1204     /* Push the bond to the bondlist */
1205     push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
1206 }
1207
1208
1209 static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr,
1210                           int atomicnumber,
1211                           int type, char *ctype, int ptype,
1212                           char *resnumberic, int cgnumber,
1213                           char *resname, char *name, real m0, real q0,
1214                           int typeB, char *ctypeB, real mB, real qB)
1215 {
1216     int           j, resind = 0, resnr;
1217     unsigned char ric;
1218     int           nr = at->nr;
1219
1220     if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1221     {
1222         gmx_fatal(FARGS, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr);
1223     }
1224
1225     j = strlen(resnumberic) - 1;
1226     if (isdigit(resnumberic[j]))
1227     {
1228         ric = ' ';
1229     }
1230     else
1231     {
1232         ric = resnumberic[j];
1233         if (j == 0 || !isdigit(resnumberic[j-1]))
1234         {
1235             gmx_fatal(FARGS, "Invalid residue number '%s' for atom %d",
1236                       resnumberic, atomnr);
1237         }
1238     }
1239     resnr = strtol(resnumberic, NULL, 10);
1240
1241     if (nr > 0)
1242     {
1243         resind = at->atom[nr-1].resind;
1244     }
1245     if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0 ||
1246         resnr != at->resinfo[resind].nr ||
1247         ric   != at->resinfo[resind].ic)
1248     {
1249         if (nr == 0)
1250         {
1251             resind = 0;
1252         }
1253         else
1254         {
1255             resind++;
1256         }
1257         at->nres = resind + 1;
1258         srenew(at->resinfo, at->nres);
1259         at->resinfo[resind].name = put_symtab(symtab, resname);
1260         at->resinfo[resind].nr   = resnr;
1261         at->resinfo[resind].ic   = ric;
1262     }
1263     else
1264     {
1265         resind = at->atom[at->nr-1].resind;
1266     }
1267
1268     /* New atom instance
1269      * get new space for arrays
1270      */
1271     srenew(at->atom, nr+1);
1272     srenew(at->atomname, nr+1);
1273     srenew(at->atomtype, nr+1);
1274     srenew(at->atomtypeB, nr+1);
1275
1276     /* fill the list */
1277     at->atom[nr].type  = type;
1278     at->atom[nr].ptype = ptype;
1279     at->atom[nr].q     = q0;
1280     at->atom[nr].m     = m0;
1281     at->atom[nr].typeB = typeB;
1282     at->atom[nr].qB    = qB;
1283     at->atom[nr].mB    = mB;
1284
1285     at->atom[nr].resind     = resind;
1286     at->atom[nr].atomnumber = atomicnumber;
1287     at->atomname[nr]        = put_symtab(symtab, name);
1288     at->atomtype[nr]        = put_symtab(symtab, ctype);
1289     at->atomtypeB[nr]       = put_symtab(symtab, ctypeB);
1290     at->nr++;
1291 }
1292
1293 void push_cg(t_block *block, int *lastindex, int index, int a)
1294 {
1295     if (debug)
1296     {
1297         fprintf (debug, "Index %d, Atom %d\n", index, a);
1298     }
1299
1300     if (((block->nr) && (*lastindex != index)) || (!block->nr))
1301     {
1302         /* add a new block */
1303         block->nr++;
1304         srenew(block->index, block->nr+1);
1305     }
1306     block->index[block->nr] = a + 1;
1307     *lastindex              = index;
1308 }
1309
1310 void push_atom(t_symtab *symtab, t_block *cgs,
1311                t_atoms *at, gpp_atomtype_t atype, char *line, int *lastcg,
1312                warninp_t wi)
1313 {
1314     int           nr, ptype;
1315     int           cgnumber, atomnr, type, typeB, nscan;
1316     char          id[STRLEN], ctype[STRLEN], ctypeB[STRLEN],
1317                   resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN];
1318     double        m, q, mb, qb;
1319     real          m0, q0, mB, qB;
1320
1321     /* Make a shortcut for writing in this molecule  */
1322     nr = at->nr;
1323
1324     /* Fixed parameters */
1325     if (sscanf(line, "%s%s%s%s%s%d",
1326                id, ctype, resnumberic, resname, name, &cgnumber) != 6)
1327     {
1328         too_few(wi);
1329         return;
1330     }
1331     sscanf(id, "%d", &atomnr);
1332     if ((type  = get_atomtype_type(ctype, atype)) == NOTSET)
1333     {
1334         gmx_fatal(FARGS, "Atomtype %s not found", ctype);
1335     }
1336     ptype = get_atomtype_ptype(type, atype);
1337
1338     /* Set default from type */
1339     q0    = get_atomtype_qA(type, atype);
1340     m0    = get_atomtype_massA(type, atype);
1341     typeB = type;
1342     qB    = q0;
1343     mB    = m0;
1344
1345     /* Optional parameters */
1346     nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1347                    &q, &m, ctypeB, &qb, &mb, check);
1348
1349     /* Nasty switch that falls thru all the way down! */
1350     if (nscan > 0)
1351     {
1352         q0 = qB = q;
1353         if (nscan > 1)
1354         {
1355             m0 = mB = m;
1356             if (nscan > 2)
1357             {
1358                 if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET)
1359                 {
1360                     gmx_fatal(FARGS, "Atomtype %s not found", ctypeB);
1361                 }
1362                 qB = get_atomtype_qA(typeB, atype);
1363                 mB = get_atomtype_massA(typeB, atype);
1364                 if (nscan > 3)
1365                 {
1366                     qB = qb;
1367                     if (nscan > 4)
1368                     {
1369                         mB = mb;
1370                         if (nscan > 5)
1371                         {
1372                             warning_error(wi, "Too many parameters");
1373                         }
1374                     }
1375                 }
1376             }
1377         }
1378     }
1379     if (debug)
1380     {
1381         fprintf(debug, "mB=%g, qB=%g, typeB=%d\n", mB, qB, typeB);
1382     }
1383
1384     push_cg(cgs, lastcg, cgnumber, nr);
1385
1386     push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype),
1387                   type, ctype, ptype, resnumberic, cgnumber,
1388                   resname, name, m0, q0, typeB,
1389                   typeB == type ? ctype : ctypeB, mB, qB);
1390 }
1391
1392 void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
1393                warninp_t wi)
1394 {
1395     char       type[STRLEN];
1396     int        nrexcl, i;
1397     t_molinfo *newmol;
1398
1399     if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
1400     {
1401         warning_error(wi, "Expected a molecule type name and nrexcl");
1402     }
1403
1404     /* Test if this atomtype overwrites another */
1405     i = 0;
1406     while (i < *nmol)
1407     {
1408         if (gmx_strcasecmp(*((*mol)[i].name), type) == 0)
1409         {
1410             gmx_fatal(FARGS, "moleculetype %s is redefined", type);
1411         }
1412         i++;
1413     }
1414
1415     (*nmol)++;
1416     srenew(*mol, *nmol);
1417     newmol = &((*mol)[*nmol-1]);
1418     init_molinfo(newmol);
1419
1420     /* Fill in the values */
1421     newmol->name     = put_symtab(symtab, type);
1422     newmol->nrexcl   = nrexcl;
1423     newmol->excl_set = FALSE;
1424 }
1425
1426 static gmx_bool default_nb_params(int ftype, t_params bt[], t_atoms *at,
1427                                   t_param *p, int c_start, gmx_bool bB, gmx_bool bGenPairs)
1428 {
1429     int          i, j, ti, tj, ntype;
1430     gmx_bool     bFound;
1431     t_param     *pi    = NULL;
1432     int          nr    = bt[ftype].nr;
1433     int          nral  = NRAL(ftype);
1434     int          nrfp  = interaction_function[ftype].nrfpA;
1435     int          nrfpB = interaction_function[ftype].nrfpB;
1436
1437     if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1438     {
1439         return TRUE;
1440     }
1441
1442     bFound = FALSE;
1443     if (bGenPairs)
1444     {
1445         /* First test the generated-pair position to save
1446          * time when we have 1000*1000 entries for e.g. OPLS...
1447          */
1448         ntype = sqrt(nr);
1449         if (bB)
1450         {
1451             ti = at->atom[p->a[0]].typeB;
1452             tj = at->atom[p->a[1]].typeB;
1453         }
1454         else
1455         {
1456             ti = at->atom[p->a[0]].type;
1457             tj = at->atom[p->a[1]].type;
1458         }
1459         pi     = &(bt[ftype].param[ntype*ti+tj]);
1460         bFound = ((ti == pi->a[0]) && (tj == pi->a[1]));
1461     }
1462
1463     /* Search explicitly if we didnt find it */
1464     if (!bFound)
1465     {
1466         for (i = 0; ((i < nr) && !bFound); i++)
1467         {
1468             pi = &(bt[ftype].param[i]);
1469             if (bB)
1470             {
1471                 for (j = 0; ((j < nral) &&
1472                              (at->atom[p->a[j]].typeB == pi->a[j])); j++)
1473                 {
1474                     ;
1475                 }
1476             }
1477             else
1478             {
1479                 for (j = 0; ((j < nral) &&
1480                              (at->atom[p->a[j]].type == pi->a[j])); j++)
1481                 {
1482                     ;
1483                 }
1484             }
1485             bFound = (j == nral);
1486         }
1487     }
1488
1489     if (bFound)
1490     {
1491         if (bB)
1492         {
1493             if (nrfp+nrfpB > MAXFORCEPARAM)
1494             {
1495                 gmx_incons("Too many force parameters");
1496             }
1497             for (j = c_start; (j < nrfpB); j++)
1498             {
1499                 p->c[nrfp+j] = pi->c[j];
1500             }
1501         }
1502         else
1503         {
1504             for (j = c_start; (j < nrfp); j++)
1505             {
1506                 p->c[j] = pi->c[j];
1507             }
1508         }
1509     }
1510     else
1511     {
1512         for (j = c_start; (j < nrfp); j++)
1513         {
1514             p->c[j] = 0.0;
1515         }
1516     }
1517     return bFound;
1518 }
1519
1520 static gmx_bool default_cmap_params(int ftype, t_params bondtype[],
1521                                     t_atoms *at, gpp_atomtype_t atype,
1522                                     t_param *p, gmx_bool bB,
1523                                     int *cmap_type, int *nparam_def)
1524 {
1525     int      i, j, nparam_found;
1526     int      ct;
1527     gmx_bool bFound = FALSE;
1528
1529     nparam_found = 0;
1530     ct           = 0;
1531
1532     /* Match the current cmap angle against the list of cmap_types */
1533     for (i = 0; i < bondtype->nct && !bFound; i += 6)
1534     {
1535         if (bB)
1536         {
1537
1538         }
1539         else
1540         {
1541             if (
1542                 (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype->cmap_types[i])   &&
1543                 (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype->cmap_types[i+1]) &&
1544                 (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype->cmap_types[i+2]) &&
1545                 (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype->cmap_types[i+3]) &&
1546                 (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype->cmap_types[i+4]))
1547             {
1548                 /* Found cmap torsion */
1549                 bFound       = TRUE;
1550                 ct           = bondtype->cmap_types[i+5];
1551                 nparam_found = 1;
1552             }
1553         }
1554     }
1555
1556     /* If we did not find a matching type for this cmap torsion */
1557     if (!bFound)
1558     {
1559         gmx_fatal(FARGS, "Unknown cmap torsion between atoms %d %d %d %d %d\n",
1560                   p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1);
1561     }
1562
1563     *nparam_def = nparam_found;
1564     *cmap_type  = ct;
1565
1566     return bFound;
1567 }
1568
1569 static gmx_bool default_params(int ftype, t_params bt[],
1570                                t_atoms *at, gpp_atomtype_t atype,
1571                                t_param *p, gmx_bool bB,
1572                                t_param **param_def,
1573                                int *nparam_def)
1574 {
1575     int          i, j, nparam_found;
1576     gmx_bool     bFound, bSame;
1577     t_param     *pi    = NULL;
1578     t_param     *pj    = NULL;
1579     int          nr    = bt[ftype].nr;
1580     int          nral  = NRAL(ftype);
1581     int          nrfpA = interaction_function[ftype].nrfpA;
1582     int          nrfpB = interaction_function[ftype].nrfpB;
1583
1584     if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1585     {
1586         return TRUE;
1587     }
1588
1589
1590     /* We allow wildcards now. The first type (with or without wildcards) that
1591      * fits is used, so you should probably put the wildcarded bondtypes
1592      * at the end of each section.
1593      */
1594     bFound       = FALSE;
1595     nparam_found = 0;
1596     /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1597      * special case for this. Check for B state outside loop to speed it up.
1598      */
1599     if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1600     {
1601         if (bB)
1602         {
1603             for (i = 0; ((i < nr) && !bFound); i++)
1604             {
1605                 pi     = &(bt[ftype].param[i]);
1606                 bFound =
1607                     (
1608                         ((pi->AI == -1) || (get_atomtype_batype(at->atom[p->AI].typeB, atype) == pi->AI)) &&
1609                         ((pi->AJ == -1) || (get_atomtype_batype(at->atom[p->AJ].typeB, atype) == pi->AJ)) &&
1610                         ((pi->AK == -1) || (get_atomtype_batype(at->atom[p->AK].typeB, atype) == pi->AK)) &&
1611                         ((pi->AL == -1) || (get_atomtype_batype(at->atom[p->AL].typeB, atype) == pi->AL))
1612                     );
1613             }
1614         }
1615         else
1616         {
1617             /* State A */
1618             for (i = 0; ((i < nr) && !bFound); i++)
1619             {
1620                 pi     = &(bt[ftype].param[i]);
1621                 bFound =
1622                     (
1623                         ((pi->AI == -1) || (get_atomtype_batype(at->atom[p->AI].type, atype) == pi->AI)) &&
1624                         ((pi->AJ == -1) || (get_atomtype_batype(at->atom[p->AJ].type, atype) == pi->AJ)) &&
1625                         ((pi->AK == -1) || (get_atomtype_batype(at->atom[p->AK].type, atype) == pi->AK)) &&
1626                         ((pi->AL == -1) || (get_atomtype_batype(at->atom[p->AL].type, atype) == pi->AL))
1627                     );
1628             }
1629         }
1630         /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1631          * The rules in that case is that additional matches HAVE to be on adjacent lines!
1632          */
1633         if (bFound == TRUE)
1634         {
1635             nparam_found++;
1636             bSame = TRUE;
1637             /* Continue from current i value */
1638             for (j = i+1; j < nr && bSame; j += 2)
1639             {
1640                 pj    = &(bt[ftype].param[j]);
1641                 bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1642                 if (bSame)
1643                 {
1644                     nparam_found++;
1645                 }
1646                 /* nparam_found will be increased as long as the numbers match */
1647             }
1648         }
1649     }
1650     else   /* Not a dihedral */
1651     {
1652         for (i = 0; ((i < nr) && !bFound); i++)
1653         {
1654             pi = &(bt[ftype].param[i]);
1655             if (bB)
1656             {
1657                 for (j = 0; ((j < nral) &&
1658                              (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++)
1659                 {
1660                     ;
1661                 }
1662             }
1663             else
1664             {
1665                 for (j = 0; ((j < nral) &&
1666                              (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++)
1667                 {
1668                     ;
1669                 }
1670             }
1671             bFound = (j == nral);
1672         }
1673         if (bFound)
1674         {
1675             nparam_found = 1;
1676         }
1677     }
1678
1679     *param_def  = pi;
1680     *nparam_def = nparam_found;
1681
1682     return bFound;
1683 }
1684
1685
1686
1687 void push_bond(directive d, t_params bondtype[], t_params bond[],
1688                t_atoms *at, gpp_atomtype_t atype, char *line,
1689                gmx_bool bBonded, gmx_bool bGenPairs, real fudgeQQ,
1690                gmx_bool bZero, gmx_bool *bWarn_copy_A_B,
1691                warninp_t wi)
1692 {
1693     const char  *aaformat[MAXATOMLIST] = {
1694         "%d%d",
1695         "%d%d%d",
1696         "%d%d%d%d",
1697         "%d%d%d%d%d",
1698         "%d%d%d%d%d%d",
1699         "%d%d%d%d%d%d%d"
1700     };
1701     const char  *asformat[MAXATOMLIST] = {
1702         "%*s%*s",
1703         "%*s%*s%*s",
1704         "%*s%*s%*s%*s",
1705         "%*s%*s%*s%*s%*s",
1706         "%*s%*s%*s%*s%*s%*s",
1707         "%*s%*s%*s%*s%*s%*s%*s"
1708     };
1709     const char  *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1710     int          nr, i, j, nral, nral_fmt, nread, ftype;
1711     char         format[STRLEN];
1712     /* One force parameter more, so we can check if we read too many */
1713     double       cc[MAXFORCEPARAM+1];
1714     int          aa[MAXATOMLIST+1];
1715     t_param      param, paramB, *param_defA, *param_defB;
1716     gmx_bool     bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE;
1717     int          nparam_defA, nparam_defB;
1718     char         errbuf[256];
1719
1720     nparam_defA = nparam_defB = 0;
1721
1722     ftype = ifunc_index(d, 1);
1723     nral  = NRAL(ftype);
1724     for (j = 0; j < MAXATOMLIST; j++)
1725     {
1726         aa[j] = NOTSET;
1727     }
1728     bDef = (NRFP(ftype) > 0);
1729
1730     if (ftype == F_SETTLE)
1731     {
1732         /* SETTLE acts on 3 atoms, but the topology format only specifies
1733          * the first atom (for historical reasons).
1734          */
1735         nral_fmt = 1;
1736     }
1737     else
1738     {
1739         nral_fmt = nral;
1740     }
1741
1742     nread = sscanf(line, aaformat[nral_fmt-1],
1743                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1744
1745     if (ftype == F_SETTLE)
1746     {
1747         aa[3] = aa[1];
1748         aa[1] = aa[0] + 1;
1749         aa[2] = aa[0] + 2;
1750     }
1751
1752     if (nread < nral_fmt)
1753     {
1754         too_few(wi);
1755         return;
1756     }
1757     else if (nread > nral_fmt)
1758     {
1759         /* this is a hack to allow for virtual sites with swapped parity */
1760         bSwapParity = (aa[nral] < 0);
1761         if (bSwapParity)
1762         {
1763             aa[nral] = -aa[nral];
1764         }
1765         ftype = ifunc_index(d, aa[nral]);
1766         if (bSwapParity)
1767         {
1768             switch (ftype)
1769             {
1770                 case F_VSITE3FAD:
1771                 case F_VSITE3OUT:
1772                     break;
1773                 default:
1774                     gmx_fatal(FARGS, "Negative function types only allowed for %s and %s",
1775                               interaction_function[F_VSITE3FAD].longname,
1776                               interaction_function[F_VSITE3OUT].longname);
1777             }
1778         }
1779     }
1780
1781
1782     /* Check for double atoms and atoms out of bounds */
1783     for (i = 0; (i < nral); i++)
1784     {
1785         if (aa[i] < 1 || aa[i] > at->nr)
1786         {
1787             gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
1788                       "Atom index (%d) in %s out of bounds (1-%d).\n"
1789                       "This probably means that you have inserted topology section \"%s\"\n"
1790                       "in a part belonging to a different molecule than you intended to.\n"
1791                       "In that case move the \"%s\" section to the right molecule.",
1792                       get_warning_file(wi), get_warning_line(wi),
1793                       aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
1794         }
1795         for (j = i+1; (j < nral); j++)
1796         {
1797             if (aa[i] == aa[j])
1798             {
1799                 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1800                 warning(wi, errbuf);
1801             }
1802         }
1803     }
1804
1805     /* default force parameters  */
1806     for (j = 0; (j < MAXATOMLIST); j++)
1807     {
1808         param.a[j] = aa[j]-1;
1809     }
1810     for (j = 0; (j < MAXFORCEPARAM); j++)
1811     {
1812         param.c[j] = 0.0;
1813     }
1814
1815     /* Get force params for normal and free energy perturbation
1816      * studies, as determined by types!
1817      */
1818
1819     if (bBonded)
1820     {
1821         bFoundA = default_params(ftype, bondtype, at, atype, &param, FALSE, &param_defA, &nparam_defA);
1822         if (bFoundA)
1823         {
1824             /* Copy the A-state and B-state default parameters */
1825             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
1826             {
1827                 param.c[j] = param_defA->c[j];
1828             }
1829         }
1830         bFoundB = default_params(ftype, bondtype, at, atype, &param, TRUE, &param_defB, &nparam_defB);
1831         if (bFoundB)
1832         {
1833             /* Copy only the B-state default parameters */
1834             for (j = NRFPA(ftype); (j < NRFP(ftype)); j++)
1835             {
1836                 param.c[j] = param_defB->c[j];
1837             }
1838         }
1839     }
1840     else if (ftype == F_LJ14)
1841     {
1842         bFoundA = default_nb_params(ftype, bondtype, at, &param, 0, FALSE, bGenPairs);
1843         bFoundB = default_nb_params(ftype, bondtype, at, &param, 0, TRUE, bGenPairs);
1844     }
1845     else if (ftype == F_LJC14_Q)
1846     {
1847         param.c[0] = fudgeQQ;
1848         /* Fill in the A-state charges as default parameters */
1849         param.c[1] = at->atom[param.a[0]].q;
1850         param.c[2] = at->atom[param.a[1]].q;
1851         /* The default LJ parameters are the standard 1-4 parameters */
1852         bFoundA = default_nb_params(F_LJ14, bondtype, at, &param, 3, FALSE, bGenPairs);
1853         bFoundB = TRUE;
1854     }
1855     else if (ftype == F_LJC_PAIRS_NB)
1856     {
1857         /* Defaults are not supported here */
1858         bFoundA = FALSE;
1859         bFoundB = TRUE;
1860     }
1861     else
1862     {
1863         gmx_incons("Unknown function type in push_bond");
1864     }
1865
1866     if (nread > nral_fmt)
1867     {
1868         /* Manually specified parameters - in this case we discard multiple torsion info! */
1869
1870         strcpy(format, asformat[nral_fmt-1]);
1871         strcat(format, ccformat);
1872
1873         nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5],
1874                        &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
1875
1876         if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1877         {
1878             /* We only have to issue a warning if these atoms are perturbed! */
1879             bPert = FALSE;
1880             for (j = 0; (j < nral); j++)
1881             {
1882                 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1883             }
1884
1885             if (bPert && *bWarn_copy_A_B)
1886             {
1887                 sprintf(errbuf,
1888                         "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1889                 warning(wi, errbuf);
1890                 *bWarn_copy_A_B = FALSE;
1891             }
1892
1893             /* If only the A parameters were specified, copy them to the B state */
1894             /* The B-state parameters correspond to the first nrfpB
1895              * A-state parameters.
1896              */
1897             for (j = 0; (j < NRFPB(ftype)); j++)
1898             {
1899                 cc[nread++] = cc[j];
1900             }
1901         }
1902
1903         /* If nread was 0 or EOF, no parameters were read => use defaults.
1904          * If nread was nrfpA we copied above so nread=nrfp.
1905          * If nread was nrfp we are cool.
1906          * For F_LJC14_Q we allow supplying fudgeQQ only.
1907          * Anything else is an error!
1908          */
1909         if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1910             !(ftype == F_LJC14_Q && nread == 1))
1911         {
1912             gmx_fatal(FARGS, "Incorrect number of parameters - found %d, expected %d or %d for %s.",
1913                       nread, NRFPA(ftype), NRFP(ftype),
1914                       interaction_function[ftype].longname);
1915         }
1916
1917         for (j = 0; (j < nread); j++)
1918         {
1919             param.c[j] = cc[j];
1920         }
1921
1922         /* Check whether we have to use the defaults */
1923         if (nread == NRFP(ftype))
1924         {
1925             bDef = FALSE;
1926         }
1927     }
1928     else
1929     {
1930         nread = 0;
1931     }
1932     /* nread now holds the number of force parameters read! */
1933
1934     if (bDef)
1935     {
1936         /* Use defaults */
1937         /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1938         if (ftype == F_PDIHS)
1939         {
1940             if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB)))
1941             {
1942                 sprintf(errbuf,
1943                         "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1944                         "Please specify perturbed parameters manually for this torsion in your topology!");
1945                 warning_error(wi, errbuf);
1946             }
1947         }
1948
1949         if (nread > 0 && nread < NRFPA(ftype))
1950         {
1951             /* Issue an error, do not use defaults */
1952             sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
1953             warning_error(wi, errbuf);
1954         }
1955
1956         if (nread == 0 || nread == EOF)
1957         {
1958             if (!bFoundA)
1959             {
1960                 if (interaction_function[ftype].flags & IF_VSITE)
1961                 {
1962                     /* set them to NOTSET, will be calculated later */
1963                     for (j = 0; (j < MAXFORCEPARAM); j++)
1964                     {
1965                         param.c[j] = NOTSET;
1966                     }
1967
1968                     if (bSwapParity)
1969                     {
1970                         param.C1 = -1; /* flag to swap parity of vsite construction */
1971                     }
1972                 }
1973                 else
1974                 {
1975                     if (bZero)
1976                     {
1977                         fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
1978                                 interaction_function[ftype].longname);
1979                     }
1980                     else
1981                     {
1982                         sprintf(errbuf, "No default %s types", interaction_function[ftype].longname);
1983                         warning_error(wi, errbuf);
1984                     }
1985                 }
1986             }
1987             else
1988             {
1989                 if (bSwapParity)
1990                 {
1991                     switch (ftype)
1992                     {
1993                         case F_VSITE3FAD:
1994                             param.C0 = 360 - param.C0;
1995                             break;
1996                         case F_VSITE3OUT:
1997                             param.C2 = -param.C2;
1998                             break;
1999                     }
2000                 }
2001             }
2002             if (!bFoundB)
2003             {
2004                 /* We only have to issue a warning if these atoms are perturbed! */
2005                 bPert = FALSE;
2006                 for (j = 0; (j < nral); j++)
2007                 {
2008                     bPert = bPert || PERTURBED(at->atom[param.a[j]]);
2009                 }
2010
2011                 if (bPert)
2012                 {
2013                     sprintf(errbuf, "No default %s types for perturbed atoms, "
2014                             "using normal values", interaction_function[ftype].longname);
2015                     warning(wi, errbuf);
2016                 }
2017             }
2018         }
2019     }
2020
2021     if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
2022         && param.c[5] != param.c[2])
2023     {
2024         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2025                   "             %s multiplicity can not be perturbed %f!=%f",
2026                   get_warning_file(wi), get_warning_line(wi),
2027                   interaction_function[ftype].longname,
2028                   param.c[2], param.c[5]);
2029     }
2030
2031     if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
2032     {
2033         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2034                   "             %s table number can not be perturbed %d!=%d",
2035                   get_warning_file(wi), get_warning_line(wi),
2036                   interaction_function[ftype].longname,
2037                   (int)(param.c[0]+0.5), (int)(param.c[2]+0.5));
2038     }
2039
2040     /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2041     if (ftype == F_RBDIHS)
2042     {
2043         nr = 0;
2044         for (i = 0; i < NRFP(ftype); i++)
2045         {
2046             if (param.c[i] != 0)
2047             {
2048                 nr++;
2049             }
2050         }
2051         if (nr == 0)
2052         {
2053             return;
2054         }
2055     }
2056
2057     /* Put the values in the appropriate arrays */
2058     add_param_to_list (&bond[ftype], &param);
2059
2060     /* Push additional torsions from FF for ftype==9 if we have them.
2061      * We have already checked that the A/B states do not differ in this case,
2062      * so we do not have to double-check that again, or the vsite stuff.
2063      * In addition, those torsions cannot be automatically perturbed.
2064      */
2065     if (bDef && ftype == F_PDIHS)
2066     {
2067         for (i = 1; i < nparam_defA; i++)
2068         {
2069             /* Advance pointer! */
2070             param_defA += 2;
2071             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
2072             {
2073                 param.c[j] = param_defA->c[j];
2074             }
2075             /* And push the next term for this torsion */
2076             add_param_to_list (&bond[ftype], &param);
2077         }
2078     }
2079 }
2080
2081 void push_cmap(directive d, t_params bondtype[], t_params bond[],
2082                t_atoms *at, gpp_atomtype_t atype, char *line,
2083                gmx_bool *bWarn_copy_A_B,
2084                warninp_t wi)
2085 {
2086     const char *aaformat[MAXATOMLIST+1] =
2087     {
2088         "%d",
2089         "%d%d",
2090         "%d%d%d",
2091         "%d%d%d%d",
2092         "%d%d%d%d%d",
2093         "%d%d%d%d%d%d",
2094         "%d%d%d%d%d%d%d"
2095     };
2096
2097     int      i, j, nr, ftype, nral, nread, ncmap_params;
2098     int      cmap_type;
2099     int      aa[MAXATOMLIST+1];
2100     char     errbuf[256];
2101     gmx_bool bFound;
2102     t_param  param, paramB, *param_defA, *param_defB;
2103
2104     ftype        = ifunc_index(d, 1);
2105     nral         = NRAL(ftype);
2106     nr           = bondtype[ftype].nr;
2107     ncmap_params = 0;
2108
2109     nread = sscanf(line, aaformat[nral-1],
2110                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2111
2112     if (nread < nral)
2113     {
2114         too_few(wi);
2115         return;
2116     }
2117     else if (nread == nral)
2118     {
2119         ftype = ifunc_index(d, 1);
2120     }
2121
2122     /* Check for double atoms and atoms out of bounds */
2123     for (i = 0; i < nral; i++)
2124     {
2125         if (aa[i] < 1 || aa[i] > at->nr)
2126         {
2127             gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2128                       "Atom index (%d) in %s out of bounds (1-%d).\n"
2129                       "This probably means that you have inserted topology section \"%s\"\n"
2130                       "in a part belonging to a different molecule than you intended to.\n"
2131                       "In that case move the \"%s\" section to the right molecule.",
2132                       get_warning_file(wi), get_warning_line(wi),
2133                       aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2134         }
2135
2136         for (j = i+1; (j < nral); j++)
2137         {
2138             if (aa[i] == aa[j])
2139             {
2140                 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2141                 warning(wi, errbuf);
2142             }
2143         }
2144     }
2145
2146     /* default force parameters  */
2147     for (j = 0; (j < MAXATOMLIST); j++)
2148     {
2149         param.a[j] = aa[j]-1;
2150     }
2151     for (j = 0; (j < MAXFORCEPARAM); j++)
2152     {
2153         param.c[j] = 0.0;
2154     }
2155
2156     /* Get the cmap type for this cmap angle */
2157     bFound = default_cmap_params(ftype, bondtype, at, atype, &param, FALSE, &cmap_type, &ncmap_params);
2158
2159     /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2160     if (bFound && ncmap_params == 1)
2161     {
2162         /* Put the values in the appropriate arrays */
2163         param.c[0] = cmap_type;
2164         add_param_to_list(&bond[ftype], &param);
2165     }
2166     else
2167     {
2168         /* This is essentially the same check as in default_cmap_params() done one more time */
2169         gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2170                   param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2171     }
2172 }
2173
2174
2175
2176 void push_vsitesn(directive d, t_params bondtype[], t_params bond[],
2177                   t_atoms *at, gpp_atomtype_t atype, char *line,
2178                   warninp_t wi)
2179 {
2180     char   *ptr;
2181     int     type, ftype, j, n, ret, nj, a;
2182     int    *atc    = NULL;
2183     double *weight = NULL, weight_tot;
2184     t_param param;
2185
2186     /* default force parameters  */
2187     for (j = 0; (j < MAXATOMLIST); j++)
2188     {
2189         param.a[j] = NOTSET;
2190     }
2191     for (j = 0; (j < MAXFORCEPARAM); j++)
2192     {
2193         param.c[j] = 0.0;
2194     }
2195
2196     ptr  = line;
2197     ret  = sscanf(ptr, "%d%n", &a, &n);
2198     ptr += n;
2199     if (ret == 0)
2200     {
2201         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2202                   "             Expected an atom index in section \"%s\"",
2203                   get_warning_file(wi), get_warning_line(wi),
2204                   dir2str(d));
2205     }
2206
2207     param.a[0] = a - 1;
2208
2209     ret   = sscanf(ptr, "%d%n", &type, &n);
2210     ptr  += n;
2211     ftype = ifunc_index(d, type);
2212
2213     weight_tot = 0;
2214     nj         = 0;
2215     do
2216     {
2217         ret  = sscanf(ptr, "%d%n", &a, &n);
2218         ptr += n;
2219         if (ret > 0)
2220         {
2221             if (nj % 20 == 0)
2222             {
2223                 srenew(atc, nj+20);
2224                 srenew(weight, nj+20);
2225             }
2226             atc[nj] = a - 1;
2227             switch (type)
2228             {
2229                 case 1:
2230                     weight[nj] = 1;
2231                     break;
2232                 case 2:
2233                     /* Here we use the A-state mass as a parameter.
2234                      * Note that the B-state mass has no influence.
2235                      */
2236                     weight[nj] = at->atom[atc[nj]].m;
2237                     break;
2238                 case 3:
2239                     weight[nj] = -1;
2240                     ret        = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2241                     ptr       += n;
2242                     if (weight[nj] < 0)
2243                     {
2244                         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2245                                   "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
2246                                   get_warning_file(wi), get_warning_line(wi),
2247                                   nj+1, atc[nj]+1);
2248                     }
2249                     break;
2250                 default:
2251                     gmx_fatal(FARGS, "Unknown vsiten type %d", type);
2252             }
2253             weight_tot += weight[nj];
2254             nj++;
2255         }
2256     }
2257     while (ret > 0);
2258
2259     if (nj == 0)
2260     {
2261         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2262                   "             Expected more than one atom index in section \"%s\"",
2263                   get_warning_file(wi), get_warning_line(wi),
2264                   dir2str(d));
2265     }
2266
2267     if (weight_tot == 0)
2268     {
2269         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2270                   "             The total mass of the construting atoms is zero",
2271                   get_warning_file(wi), get_warning_line(wi));
2272     }
2273
2274     for (j = 0; j < nj; j++)
2275     {
2276         param.a[1] = atc[j];
2277         param.c[0] = nj;
2278         param.c[1] = weight[j]/weight_tot;
2279         /* Put the values in the appropriate arrays */
2280         add_param_to_list (&bond[ftype], &param);
2281     }
2282
2283     sfree(atc);
2284     sfree(weight);
2285 }
2286
2287 void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
2288               int *nrcopies,
2289               warninp_t wi)
2290 {
2291     int  i, copies;
2292     char type[STRLEN];
2293
2294     *nrcopies = 0;
2295     if (sscanf(pline, "%s%d", type, &copies) != 2)
2296     {
2297         too_few(wi);
2298         return;
2299     }
2300
2301     /* search moleculename */
2302     for (i = 0; ((i < nrmols) && gmx_strcasecmp(type, *(mols[i].name))); i++)
2303     {
2304         ;
2305     }
2306
2307     if (i < nrmols)
2308     {
2309         *nrcopies        = copies;
2310         *whichmol        = i;
2311     }
2312     else
2313     {
2314         gmx_fatal(FARGS, "No such moleculetype %s", type);
2315     }
2316 }
2317
2318 void init_block2(t_block2 *b2, int natom)
2319 {
2320     int i;
2321
2322     b2->nr = natom;
2323     snew(b2->nra, b2->nr);
2324     snew(b2->a, b2->nr);
2325     for (i = 0; (i < b2->nr); i++)
2326     {
2327         b2->a[i] = NULL;
2328     }
2329 }
2330
2331 void done_block2(t_block2 *b2)
2332 {
2333     int i;
2334
2335     if (b2->nr)
2336     {
2337         for (i = 0; (i < b2->nr); i++)
2338         {
2339             sfree(b2->a[i]);
2340         }
2341         sfree(b2->a);
2342         sfree(b2->nra);
2343         b2->nr = 0;
2344     }
2345 }
2346
2347 void push_excl(char *line, t_block2 *b2)
2348 {
2349     int  i, j;
2350     int  n;
2351     char base[STRLEN], format[STRLEN];
2352
2353     if (sscanf(line, "%d", &i) == 0)
2354     {
2355         return;
2356     }
2357
2358     if ((1 <= i) && (i <= b2->nr))
2359     {
2360         i--;
2361     }
2362     else
2363     {
2364         if (debug)
2365         {
2366             fprintf(debug, "Unbound atom %d\n", i-1);
2367         }
2368         return;
2369     }
2370     strcpy(base, "%*d");
2371     do
2372     {
2373         strcpy(format, base);
2374         strcat(format, "%d");
2375         n = sscanf(line, format, &j);
2376         if (n == 1)
2377         {
2378             if ((1 <= j) && (j <= b2->nr))
2379             {
2380                 j--;
2381                 srenew(b2->a[i], ++(b2->nra[i]));
2382                 b2->a[i][b2->nra[i]-1] = j;
2383                 /* also add the reverse exclusion! */
2384                 srenew(b2->a[j], ++(b2->nra[j]));
2385                 b2->a[j][b2->nra[j]-1] = i;
2386                 strcat(base, "%*d");
2387             }
2388             else
2389             {
2390                 gmx_fatal(FARGS, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
2391             }
2392         }
2393     }
2394     while (n == 1);
2395 }
2396
2397 void b_to_b2(t_blocka *b, t_block2 *b2)
2398 {
2399     int     i;
2400     atom_id j, a;
2401
2402     for (i = 0; (i < b->nr); i++)
2403     {
2404         for (j = b->index[i]; (j < b->index[i+1]); j++)
2405         {
2406             a = b->a[j];
2407             srenew(b2->a[i], ++b2->nra[i]);
2408             b2->a[i][b2->nra[i]-1] = a;
2409         }
2410     }
2411 }
2412
2413 void b2_to_b(t_block2 *b2, t_blocka *b)
2414 {
2415     int     i, nra;
2416     atom_id j;
2417
2418     nra = 0;
2419     for (i = 0; (i < b2->nr); i++)
2420     {
2421         b->index[i] = nra;
2422         for (j = 0; (j < b2->nra[i]); j++)
2423         {
2424             b->a[nra+j] = b2->a[i][j];
2425         }
2426         nra += b2->nra[i];
2427     }
2428     /* terminate list */
2429     b->index[i] = nra;
2430 }
2431
2432 static int icomp(const void *v1, const void *v2)
2433 {
2434     return (*((atom_id *) v1))-(*((atom_id *) v2));
2435 }
2436
2437 void merge_excl(t_blocka *excl, t_block2 *b2)
2438 {
2439     int     i, k;
2440     atom_id j;
2441     int     nra;
2442
2443     if (!b2->nr)
2444     {
2445         return;
2446     }
2447     else if (b2->nr != excl->nr)
2448     {
2449         gmx_fatal(FARGS, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2450                   b2->nr, excl->nr);
2451     }
2452     else if (debug)
2453     {
2454         fprintf(debug, "Entering merge_excl\n");
2455     }
2456
2457     /* First copy all entries from excl to b2 */
2458     b_to_b2(excl, b2);
2459
2460     /* Count and sort the exclusions */
2461     nra = 0;
2462     for (i = 0; (i < b2->nr); i++)
2463     {
2464         if (b2->nra[i] > 0)
2465         {
2466             /* remove double entries */
2467             qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp);
2468             k = 1;
2469             for (j = 1; (j < b2->nra[i]); j++)
2470             {
2471                 if (b2->a[i][j] != b2->a[i][k-1])
2472                 {
2473                     b2->a[i][k] = b2->a[i][j];
2474                     k++;
2475                 }
2476             }
2477             b2->nra[i] = k;
2478             nra       += b2->nra[i];
2479         }
2480     }
2481     excl->nra = nra;
2482     srenew(excl->a, excl->nra);
2483
2484     b2_to_b(b2, excl);
2485 }
2486
2487 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
2488                            t_nbparam ***nbparam, t_nbparam ***pair)
2489 {
2490     t_atom  atom;
2491     t_param param;
2492     int     i, nr;
2493
2494     /* Define an atom type with all parameters set to zero (no interactions) */
2495     atom.q = 0.0;
2496     atom.m = 0.0;
2497     /* Type for decoupled atoms could be anything,
2498      * this should be changed automatically later when required.
2499      */
2500     atom.ptype = eptAtom;
2501     for (i = 0; (i < MAXFORCEPARAM); i++)
2502     {
2503         param.c[i] = 0.0;
2504     }
2505
2506     nr = add_atomtype(at, symtab, &atom, "decoupled", &param, -1, 0.0, 0.0, 0.0, 0, 0, 0);
2507
2508     /* Add space in the non-bonded parameters matrix */
2509     realloc_nb_params(at, nbparam, pair);
2510
2511     return nr;
2512 }
2513
2514 static void convert_pairs_to_pairsQ(t_params *plist,
2515                                     real fudgeQQ, t_atoms *atoms)
2516 {
2517     t_param *param;
2518     int      i;
2519     real     v, w;
2520
2521     /* Copy the pair list to the pairQ list */
2522     plist[F_LJC14_Q] = plist[F_LJ14];
2523     /* Empty the pair list */
2524     plist[F_LJ14].nr    = 0;
2525     plist[F_LJ14].param = NULL;
2526     param               = plist[F_LJC14_Q].param;
2527     for (i = 0; i < plist[F_LJC14_Q].nr; i++)
2528     {
2529         v             = param[i].c[0];
2530         w             = param[i].c[1];
2531         param[i].c[0] = fudgeQQ;
2532         param[i].c[1] = atoms->atom[param[i].a[0]].q;
2533         param[i].c[2] = atoms->atom[param[i].a[1]].q;
2534         param[i].c[3] = v;
2535         param[i].c[4] = w;
2536     }
2537 }
2538
2539 static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp)
2540 {
2541     int       n, ntype, i, j, k;
2542     t_atom   *atom;
2543     t_blocka *excl;
2544     gmx_bool  bExcl;
2545     t_param   param;
2546
2547     n    = mol->atoms.nr;
2548     atom = mol->atoms.atom;
2549
2550     ntype = sqrt(nbp->nr);
2551
2552     for (i = 0; i < MAXATOMLIST; i++)
2553     {
2554         param.a[i] = NOTSET;
2555     }
2556     for (i = 0; i < MAXFORCEPARAM; i++)
2557     {
2558         param.c[i] = NOTSET;
2559     }
2560
2561     /* Add a pair interaction for all non-excluded atom pairs */
2562     excl = &mol->excls;
2563     for (i = 0; i < n; i++)
2564     {
2565         for (j = i+1; j < n; j++)
2566         {
2567             bExcl = FALSE;
2568             for (k = excl->index[i]; k < excl->index[i+1]; k++)
2569             {
2570                 if (excl->a[k] == j)
2571                 {
2572                     bExcl = TRUE;
2573                 }
2574             }
2575             if (!bExcl)
2576             {
2577                 if (nb_funct != F_LJ)
2578                 {
2579                     gmx_fatal(FARGS, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2580                 }
2581                 param.a[0] = i;
2582                 param.a[1] = j;
2583                 param.c[0] = atom[i].q;
2584                 param.c[1] = atom[j].q;
2585                 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2586                 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2587                 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], &param);
2588             }
2589         }
2590     }
2591 }
2592
2593 static void set_excl_all(t_blocka *excl)
2594 {
2595     int nat, i, j, k;
2596
2597     /* Get rid of the current exclusions and exclude all atom pairs */
2598     nat       = excl->nr;
2599     excl->nra = nat*nat;
2600     srenew(excl->a, excl->nra);
2601     k = 0;
2602     for (i = 0; i < nat; i++)
2603     {
2604         excl->index[i] = k;
2605         for (j = 0; j < nat; j++)
2606         {
2607             excl->a[k++] = j;
2608         }
2609     }
2610     excl->index[nat] = k;
2611 }
2612
2613 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2614                            int couple_lam0, int couple_lam1)
2615 {
2616     int i;
2617
2618     for (i = 0; i < atoms->nr; i++)
2619     {
2620         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2621         {
2622             atoms->atom[i].q     = 0.0;
2623         }
2624         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2625         {
2626             atoms->atom[i].type  = atomtype_decouple;
2627         }
2628         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2629         {
2630             atoms->atom[i].qB    = 0.0;
2631         }
2632         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2633         {
2634             atoms->atom[i].typeB = atomtype_decouple;
2635         }
2636     }
2637 }
2638
2639 void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
2640                             int couple_lam0, int couple_lam1,
2641                             gmx_bool bCoupleIntra, int nb_funct, t_params *nbp)
2642 {
2643     convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2644     if (!bCoupleIntra)
2645     {
2646         generate_LJCpairsNB(mol, nb_funct, nbp);
2647         set_excl_all(&mol->excls);
2648     }
2649     decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1);
2650 }