src/gromacs/gmxpreprocess/toppush.c

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