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