src/kernel/toppush.c

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