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