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