src/gromacs/gmxpreprocess/toppush.c

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