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,
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,
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(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                warninp_t wi)
2084 {
2085     const char *aaformat[MAXATOMLIST+1] =
2086     {
2087         "%d",
2088         "%d%d",
2089         "%d%d%d",
2090         "%d%d%d%d",
2091         "%d%d%d%d%d",
2092         "%d%d%d%d%d%d",
2093         "%d%d%d%d%d%d%d"
2094     };
2095
2096     int      i, j, nr, ftype, nral, nread, ncmap_params;
2097     int      cmap_type;
2098     int      aa[MAXATOMLIST+1];
2099     char     errbuf[256];
2100     gmx_bool bFound;
2101     t_param  param, paramB, *param_defA, *param_defB;
2102
2103     ftype        = ifunc_index(d, 1);
2104     nral         = NRAL(ftype);
2105     nr           = bondtype[ftype].nr;
2106     ncmap_params = 0;
2107
2108     nread = sscanf(line, aaformat[nral-1],
2109                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2110
2111     if (nread < nral)
2112     {
2113         too_few(wi);
2114         return;
2115     }
2116     else if (nread == nral)
2117     {
2118         ftype = ifunc_index(d, 1);
2119     }
2120
2121     /* Check for double atoms and atoms out of bounds */
2122     for (i = 0; i < nral; i++)
2123     {
2124         if (aa[i] < 1 || aa[i] > at->nr)
2125         {
2126             gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2127                       "Atom index (%d) in %s out of bounds (1-%d).\n"
2128                       "This probably means that you have inserted topology section \"%s\"\n"
2129                       "in a part belonging to a different molecule than you intended to.\n"
2130                       "In that case move the \"%s\" section to the right molecule.",
2131                       get_warning_file(wi), get_warning_line(wi),
2132                       aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2133         }
2134
2135         for (j = i+1; (j < nral); j++)
2136         {
2137             if (aa[i] == aa[j])
2138             {
2139                 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2140                 warning(wi, errbuf);
2141             }
2142         }
2143     }
2144
2145     /* default force parameters  */
2146     for (j = 0; (j < MAXATOMLIST); j++)
2147     {
2148         param.a[j] = aa[j]-1;
2149     }
2150     for (j = 0; (j < MAXFORCEPARAM); j++)
2151     {
2152         param.c[j] = 0.0;
2153     }
2154
2155     /* Get the cmap type for this cmap angle */
2156     bFound = default_cmap_params(bondtype, at, atype, &param, FALSE, &cmap_type, &ncmap_params);
2157
2158     /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2159     if (bFound && ncmap_params == 1)
2160     {
2161         /* Put the values in the appropriate arrays */
2162         param.c[0] = cmap_type;
2163         add_param_to_list(&bond[ftype], &param);
2164     }
2165     else
2166     {
2167         /* This is essentially the same check as in default_cmap_params() done one more time */
2168         gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2169                   param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2170     }
2171 }
2172
2173
2174
2175 void push_vsitesn(directive d, t_params bond[],
2176                   t_atoms *at, char *line,
2177                   warninp_t wi)
2178 {
2179     char   *ptr;
2180     int     type, ftype, j, n, ret, nj, a;
2181     int    *atc    = NULL;
2182     double *weight = NULL, weight_tot;
2183     t_param param;
2184
2185     /* default force parameters  */
2186     for (j = 0; (j < MAXATOMLIST); j++)
2187     {
2188         param.a[j] = NOTSET;
2189     }
2190     for (j = 0; (j < MAXFORCEPARAM); j++)
2191     {
2192         param.c[j] = 0.0;
2193     }
2194
2195     ptr  = line;
2196     ret  = sscanf(ptr, "%d%n", &a, &n);
2197     ptr += n;
2198     if (ret == 0)
2199     {
2200         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2201                   "             Expected an atom index in section \"%s\"",
2202                   get_warning_file(wi), get_warning_line(wi),
2203                   dir2str(d));
2204     }
2205
2206     param.a[0] = a - 1;
2207
2208     ret   = sscanf(ptr, "%d%n", &type, &n);
2209     ptr  += n;
2210     ftype = ifunc_index(d, type);
2211
2212     weight_tot = 0;
2213     nj         = 0;
2214     do
2215     {
2216         ret  = sscanf(ptr, "%d%n", &a, &n);
2217         ptr += n;
2218         if (ret > 0)
2219         {
2220             if (nj % 20 == 0)
2221             {
2222                 srenew(atc, nj+20);
2223                 srenew(weight, nj+20);
2224             }
2225             atc[nj] = a - 1;
2226             switch (type)
2227             {
2228                 case 1:
2229                     weight[nj] = 1;
2230                     break;
2231                 case 2:
2232                     /* Here we use the A-state mass as a parameter.
2233                      * Note that the B-state mass has no influence.
2234                      */
2235                     weight[nj] = at->atom[atc[nj]].m;
2236                     break;
2237                 case 3:
2238                     weight[nj] = -1;
2239                     ret        = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2240                     ptr       += n;
2241                     if (weight[nj] < 0)
2242                     {
2243                         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2244                                   "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
2245                                   get_warning_file(wi), get_warning_line(wi),
2246                                   nj+1, atc[nj]+1);
2247                     }
2248                     break;
2249                 default:
2250                     gmx_fatal(FARGS, "Unknown vsiten type %d", type);
2251             }
2252             weight_tot += weight[nj];
2253             nj++;
2254         }
2255     }
2256     while (ret > 0);
2257
2258     if (nj == 0)
2259     {
2260         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2261                   "             Expected more than one atom index in section \"%s\"",
2262                   get_warning_file(wi), get_warning_line(wi),
2263                   dir2str(d));
2264     }
2265
2266     if (weight_tot == 0)
2267     {
2268         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2269                   "             The total mass of the construting atoms is zero",
2270                   get_warning_file(wi), get_warning_line(wi));
2271     }
2272
2273     for (j = 0; j < nj; j++)
2274     {
2275         param.a[1] = atc[j];
2276         param.c[0] = nj;
2277         param.c[1] = weight[j]/weight_tot;
2278         /* Put the values in the appropriate arrays */
2279         add_param_to_list (&bond[ftype], &param);
2280     }
2281
2282     sfree(atc);
2283     sfree(weight);
2284 }
2285
2286 void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
2287               int *nrcopies,
2288               warninp_t wi)
2289 {
2290     int  i, copies;
2291     char type[STRLEN];
2292
2293     *nrcopies = 0;
2294     if (sscanf(pline, "%s%d", type, &copies) != 2)
2295     {
2296         too_few(wi);
2297         return;
2298     }
2299
2300     /* search moleculename */
2301     for (i = 0; ((i < nrmols) && gmx_strcasecmp(type, *(mols[i].name))); i++)
2302     {
2303         ;
2304     }
2305
2306     if (i < nrmols)
2307     {
2308         *nrcopies        = copies;
2309         *whichmol        = i;
2310     }
2311     else
2312     {
2313         gmx_fatal(FARGS, "No such moleculetype %s", type);
2314     }
2315 }
2316
2317 void init_block2(t_block2 *b2, int natom)
2318 {
2319     int i;
2320
2321     b2->nr = natom;
2322     snew(b2->nra, b2->nr);
2323     snew(b2->a, b2->nr);
2324     for (i = 0; (i < b2->nr); i++)
2325     {
2326         b2->a[i] = NULL;
2327     }
2328 }
2329
2330 void done_block2(t_block2 *b2)
2331 {
2332     int i;
2333
2334     if (b2->nr)
2335     {
2336         for (i = 0; (i < b2->nr); i++)
2337         {
2338             sfree(b2->a[i]);
2339         }
2340         sfree(b2->a);
2341         sfree(b2->nra);
2342         b2->nr = 0;
2343     }
2344 }
2345
2346 void push_excl(char *line, t_block2 *b2)
2347 {
2348     int  i, j;
2349     int  n;
2350     char base[STRLEN], format[STRLEN];
2351
2352     if (sscanf(line, "%d", &i) == 0)
2353     {
2354         return;
2355     }
2356
2357     if ((1 <= i) && (i <= b2->nr))
2358     {
2359         i--;
2360     }
2361     else
2362     {
2363         if (debug)
2364         {
2365             fprintf(debug, "Unbound atom %d\n", i-1);
2366         }
2367         return;
2368     }
2369     strcpy(base, "%*d");
2370     do
2371     {
2372         strcpy(format, base);
2373         strcat(format, "%d");
2374         n = sscanf(line, format, &j);
2375         if (n == 1)
2376         {
2377             if ((1 <= j) && (j <= b2->nr))
2378             {
2379                 j--;
2380                 srenew(b2->a[i], ++(b2->nra[i]));
2381                 b2->a[i][b2->nra[i]-1] = j;
2382                 /* also add the reverse exclusion! */
2383                 srenew(b2->a[j], ++(b2->nra[j]));
2384                 b2->a[j][b2->nra[j]-1] = i;
2385                 strcat(base, "%*d");
2386             }
2387             else
2388             {
2389                 gmx_fatal(FARGS, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
2390             }
2391         }
2392     }
2393     while (n == 1);
2394 }
2395
2396 void b_to_b2(t_blocka *b, t_block2 *b2)
2397 {
2398     int     i;
2399     atom_id j, a;
2400
2401     for (i = 0; (i < b->nr); i++)
2402     {
2403         for (j = b->index[i]; (j < b->index[i+1]); j++)
2404         {
2405             a = b->a[j];
2406             srenew(b2->a[i], ++b2->nra[i]);
2407             b2->a[i][b2->nra[i]-1] = a;
2408         }
2409     }
2410 }
2411
2412 void b2_to_b(t_block2 *b2, t_blocka *b)
2413 {
2414     int     i, nra;
2415     atom_id j;
2416
2417     nra = 0;
2418     for (i = 0; (i < b2->nr); i++)
2419     {
2420         b->index[i] = nra;
2421         for (j = 0; (j < b2->nra[i]); j++)
2422         {
2423             b->a[nra+j] = b2->a[i][j];
2424         }
2425         nra += b2->nra[i];
2426     }
2427     /* terminate list */
2428     b->index[i] = nra;
2429 }
2430
2431 static int icomp(const void *v1, const void *v2)
2432 {
2433     return (*((atom_id *) v1))-(*((atom_id *) v2));
2434 }
2435
2436 void merge_excl(t_blocka *excl, t_block2 *b2)
2437 {
2438     int     i, k;
2439     atom_id j;
2440     int     nra;
2441
2442     if (!b2->nr)
2443     {
2444         return;
2445     }
2446     else if (b2->nr != excl->nr)
2447     {
2448         gmx_fatal(FARGS, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2449                   b2->nr, excl->nr);
2450     }
2451     else if (debug)
2452     {
2453         fprintf(debug, "Entering merge_excl\n");
2454     }
2455
2456     /* First copy all entries from excl to b2 */
2457     b_to_b2(excl, b2);
2458
2459     /* Count and sort the exclusions */
2460     nra = 0;
2461     for (i = 0; (i < b2->nr); i++)
2462     {
2463         if (b2->nra[i] > 0)
2464         {
2465             /* remove double entries */
2466             qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp);
2467             k = 1;
2468             for (j = 1; (j < b2->nra[i]); j++)
2469             {
2470                 if (b2->a[i][j] != b2->a[i][k-1])
2471                 {
2472                     b2->a[i][k] = b2->a[i][j];
2473                     k++;
2474                 }
2475             }
2476             b2->nra[i] = k;
2477             nra       += b2->nra[i];
2478         }
2479     }
2480     excl->nra = nra;
2481     srenew(excl->a, excl->nra);
2482
2483     b2_to_b(b2, excl);
2484 }
2485
2486 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
2487                            t_nbparam ***nbparam, t_nbparam ***pair)
2488 {
2489     t_atom  atom;
2490     t_param param;
2491     int     i, nr;
2492
2493     /* Define an atom type with all parameters set to zero (no interactions) */
2494     atom.q = 0.0;
2495     atom.m = 0.0;
2496     /* Type for decoupled atoms could be anything,
2497      * this should be changed automatically later when required.
2498      */
2499     atom.ptype = eptAtom;
2500     for (i = 0; (i < MAXFORCEPARAM); i++)
2501     {
2502         param.c[i] = 0.0;
2503     }
2504
2505     nr = add_atomtype(at, symtab, &atom, "decoupled", &param, -1, 0.0, 0.0, 0.0, 0, 0, 0);
2506
2507     /* Add space in the non-bonded parameters matrix */
2508     realloc_nb_params(at, nbparam, pair);
2509
2510     return nr;
2511 }
2512
2513 static void convert_pairs_to_pairsQ(t_params *plist,
2514                                     real fudgeQQ, t_atoms *atoms)
2515 {
2516     t_param *param;
2517     int      i;
2518     real     v, w;
2519
2520     /* Copy the pair list to the pairQ list */
2521     plist[F_LJC14_Q] = plist[F_LJ14];
2522     /* Empty the pair list */
2523     plist[F_LJ14].nr    = 0;
2524     plist[F_LJ14].param = NULL;
2525     param               = plist[F_LJC14_Q].param;
2526     for (i = 0; i < plist[F_LJC14_Q].nr; i++)
2527     {
2528         v             = param[i].c[0];
2529         w             = param[i].c[1];
2530         param[i].c[0] = fudgeQQ;
2531         param[i].c[1] = atoms->atom[param[i].a[0]].q;
2532         param[i].c[2] = atoms->atom[param[i].a[1]].q;
2533         param[i].c[3] = v;
2534         param[i].c[4] = w;
2535     }
2536 }
2537
2538 static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp)
2539 {
2540     int       n, ntype, i, j, k;
2541     t_atom   *atom;
2542     t_blocka *excl;
2543     gmx_bool  bExcl;
2544     t_param   param;
2545
2546     n    = mol->atoms.nr;
2547     atom = mol->atoms.atom;
2548
2549     ntype = sqrt(nbp->nr);
2550
2551     for (i = 0; i < MAXATOMLIST; i++)
2552     {
2553         param.a[i] = NOTSET;
2554     }
2555     for (i = 0; i < MAXFORCEPARAM; i++)
2556     {
2557         param.c[i] = NOTSET;
2558     }
2559
2560     /* Add a pair interaction for all non-excluded atom pairs */
2561     excl = &mol->excls;
2562     for (i = 0; i < n; i++)
2563     {
2564         for (j = i+1; j < n; j++)
2565         {
2566             bExcl = FALSE;
2567             for (k = excl->index[i]; k < excl->index[i+1]; k++)
2568             {
2569                 if (excl->a[k] == j)
2570                 {
2571                     bExcl = TRUE;
2572                 }
2573             }
2574             if (!bExcl)
2575             {
2576                 if (nb_funct != F_LJ)
2577                 {
2578                     gmx_fatal(FARGS, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2579                 }
2580                 param.a[0] = i;
2581                 param.a[1] = j;
2582                 param.c[0] = atom[i].q;
2583                 param.c[1] = atom[j].q;
2584                 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2585                 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2586                 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], &param);
2587             }
2588         }
2589     }
2590 }
2591
2592 static void set_excl_all(t_blocka *excl)
2593 {
2594     int nat, i, j, k;
2595
2596     /* Get rid of the current exclusions and exclude all atom pairs */
2597     nat       = excl->nr;
2598     excl->nra = nat*nat;
2599     srenew(excl->a, excl->nra);
2600     k = 0;
2601     for (i = 0; i < nat; i++)
2602     {
2603         excl->index[i] = k;
2604         for (j = 0; j < nat; j++)
2605         {
2606             excl->a[k++] = j;
2607         }
2608     }
2609     excl->index[nat] = k;
2610 }
2611
2612 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2613                            int couple_lam0, int couple_lam1)
2614 {
2615     int i;
2616
2617     for (i = 0; i < atoms->nr; i++)
2618     {
2619         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2620         {
2621             atoms->atom[i].q     = 0.0;
2622         }
2623         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2624         {
2625             atoms->atom[i].type  = atomtype_decouple;
2626         }
2627         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2628         {
2629             atoms->atom[i].qB    = 0.0;
2630         }
2631         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2632         {
2633             atoms->atom[i].typeB = atomtype_decouple;
2634         }
2635     }
2636 }
2637
2638 void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
2639                             int couple_lam0, int couple_lam1,
2640                             gmx_bool bCoupleIntra, int nb_funct, t_params *nbp)
2641 {
2642     convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2643     if (!bCoupleIntra)
2644     {
2645         generate_LJCpairsNB(mol, nb_funct, nbp);
2646         set_excl_all(&mol->excls);
2647     }
2648     decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1);
2649 }