src/gromacs/gmxpreprocess/toppush.c

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