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