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