src/gromacs/gmxpreprocess/toppush.cpp

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