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, InteractionTypeParameters *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 InteractionTypeParameters. */
 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, InteractionTypeParameters *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(InteractionTypeParameters     *       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, gmx::ArrayRef<InteractionTypeParameters> 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, gmx::ArrayRef<InteractionTypeParameters> 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, gmx::ArrayRef<InteractionTypeParameters> 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     /* Read in CMAP parameters */
1159     sl = 0;
1160     for (i = 0; i < ncmap; i++)
1161     {
1162         while (isspace(*(line+start+sl)))
1163         {
1164             sl++;
1165         }
1166         nn  = sscanf(line+start+sl, " %s ", s);
1167         sl += strlen(s);
1168         bt[F_CMAP].cmap.emplace_back(strtod(s, nullptr));
1169
1170         if (nn == 1)
1171         {
1172             read_cmap++;
1173         }
1174         else
1175         {
1176             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]);
1177             warning_error(wi, message);
1178         }
1179
1180     }
1181
1182     /* Check do that we got the number of parameters we expected */
1183     if (read_cmap == nrfpA)
1184     {
1185         for (i = 0; i < ncmap; i++)
1186         {
1187             bt[F_CMAP].cmap.emplace_back(bt[F_CMAP].cmap[i]);
1188         }
1189     }
1190     else
1191     {
1192         if (read_cmap < nrfpA)
1193         {
1194             warning_error(wi, "Not enough cmap parameters");
1195         }
1196         else if (read_cmap > nrfpA && read_cmap < nrfp)
1197         {
1198             warning_error(wi, "Too many cmap parameters or not enough parameters for topology B");
1199         }
1200         else if (read_cmap > nrfp)
1201         {
1202             warning_error(wi, "Too many cmap parameters");
1203         }
1204     }
1205
1206
1207     /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1208      * so we can safely assign them each time
1209      */
1210     bt[F_CMAP].cmakeGridSpacing = nxcmap; /* Or nycmap, they need to be equal */
1211     bt[F_CMAP].cmapAngles++;              /* Since we are incrementing here, we need to subtract later, see (*****) */
1212     nct                     = (nral+1) * bt[F_CMAP].cmapAngles;
1213
1214     for (i = 0; (i < nral); i++)
1215     {
1216         if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1217         {
1218             auto message = gmx::formatString("Unknown atomtype %s\n", alc[i]);
1219             warning_error(wi, message);
1220         }
1221         else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1222         {
1223             auto message = gmx::formatString("Unknown bond_atomtype %s\n", alc[i]);
1224             warning_error(wi, message);
1225         }
1226
1227         /* Assign a grid number to each cmap_type */
1228         bt[F_CMAP].cmapAtomTypes.emplace_back(get_bond_atomtype_type(alc[i], bat));
1229     }
1230
1231     /* Assign a type number to this cmap */
1232     bt[F_CMAP].cmapAtomTypes.emplace_back(bt[F_CMAP].cmapAngles-1); /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1233
1234     /* Check for the correct number of atoms (again) */
1235     if (bt[F_CMAP].nct() != nct)
1236     {
1237         auto message = gmx::formatString("Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].cmapAngles);
1238         warning_error(wi, message);
1239     }
1240
1241     /* Is this correct?? */
1242     for (i = 0; i < MAXFORCEPARAM; i++)
1243     {
1244         p.c[i] = NOTSET;
1245     }
1246
1247     /* Push the bond to the bondlist */
1248     push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
1249 }
1250
1251
1252 static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr,
1253                           int atomicnumber,
1254                           int type, char *ctype, int ptype,
1255                           char *resnumberic,
1256                           char *resname, char *name, real m0, real q0,
1257                           int typeB, char *ctypeB, real mB, real qB,
1258                           warninp *wi)
1259 {
1260     int           j, resind = 0, resnr;
1261     unsigned char ric;
1262     int           nr = at->nr;
1263
1264     if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1265     {
1266         auto message = gmx::formatString
1267                 ("Atoms in the .top are not numbered consecutively from 1 (rather, "
1268                 "atomnr = %d, while at->nr = %d)", atomnr, at->nr);
1269         warning_error_and_exit(wi, message, FARGS);
1270     }
1271
1272     j = strlen(resnumberic) - 1;
1273     if (isdigit(resnumberic[j]))
1274     {
1275         ric = ' ';
1276     }
1277     else
1278     {
1279         ric = resnumberic[j];
1280         if (j == 0 || !isdigit(resnumberic[j-1]))
1281         {
1282             auto message = gmx::formatString("Invalid residue number '%s' for atom %d",
1283                                              resnumberic, atomnr);
1284             warning_error_and_exit(wi, message, FARGS);
1285         }
1286     }
1287     resnr = strtol(resnumberic, nullptr, 10);
1288
1289     if (nr > 0)
1290     {
1291         resind = at->atom[nr-1].resind;
1292     }
1293     if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0 ||
1294         resnr != at->resinfo[resind].nr ||
1295         ric   != at->resinfo[resind].ic)
1296     {
1297         if (nr == 0)
1298         {
1299             resind = 0;
1300         }
1301         else
1302         {
1303             resind++;
1304         }
1305         at->nres = resind + 1;
1306         srenew(at->resinfo, at->nres);
1307         at->resinfo[resind].name = put_symtab(symtab, resname);
1308         at->resinfo[resind].nr   = resnr;
1309         at->resinfo[resind].ic   = ric;
1310     }
1311     else
1312     {
1313         resind = at->atom[at->nr-1].resind;
1314     }
1315
1316     /* New atom instance
1317      * get new space for arrays
1318      */
1319     srenew(at->atom, nr+1);
1320     srenew(at->atomname, nr+1);
1321     srenew(at->atomtype, nr+1);
1322     srenew(at->atomtypeB, nr+1);
1323
1324     /* fill the list */
1325     at->atom[nr].type  = type;
1326     at->atom[nr].ptype = ptype;
1327     at->atom[nr].q     = q0;
1328     at->atom[nr].m     = m0;
1329     at->atom[nr].typeB = typeB;
1330     at->atom[nr].qB    = qB;
1331     at->atom[nr].mB    = mB;
1332
1333     at->atom[nr].resind     = resind;
1334     at->atom[nr].atomnumber = atomicnumber;
1335     at->atomname[nr]        = put_symtab(symtab, name);
1336     at->atomtype[nr]        = put_symtab(symtab, ctype);
1337     at->atomtypeB[nr]       = put_symtab(symtab, ctypeB);
1338     at->nr++;
1339 }
1340
1341 static void push_cg(t_block *block, int *lastindex, int index, int a)
1342 {
1343     if (((block->nr) && (*lastindex != index)) || (!block->nr))
1344     {
1345         /* add a new block */
1346         block->nr++;
1347         srenew(block->index, block->nr+1);
1348     }
1349     block->index[block->nr] = a + 1;
1350     *lastindex              = index;
1351 }
1352
1353 void push_atom(t_symtab *symtab, t_block *cgs,
1354                t_atoms *at, gpp_atomtype *atype, char *line, int *lastcg,
1355                warninp *wi)
1356 {
1357     int           nr, ptype;
1358     int           cgnumber, atomnr, type, typeB, nscan;
1359     char          id[STRLEN], ctype[STRLEN], ctypeB[STRLEN],
1360                   resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN];
1361     double        m, q, mb, qb;
1362     real          m0, q0, mB, qB;
1363
1364     /* Make a shortcut for writing in this molecule  */
1365     nr = at->nr;
1366
1367     /* Fixed parameters */
1368     if (sscanf(line, "%s%s%s%s%s%d",
1369                id, ctype, resnumberic, resname, name, &cgnumber) != 6)
1370     {
1371         too_few(wi);
1372         return;
1373     }
1374     sscanf(id, "%d", &atomnr);
1375     if ((type  = get_atomtype_type(ctype, atype)) == NOTSET)
1376     {
1377         auto message = gmx::formatString("Atomtype %s not found", ctype);
1378         warning_error_and_exit(wi, message, FARGS);
1379     }
1380     ptype = get_atomtype_ptype(type, atype);
1381
1382     /* Set default from type */
1383     q0    = get_atomtype_qA(type, atype);
1384     m0    = get_atomtype_massA(type, atype);
1385     typeB = type;
1386     qB    = q0;
1387     mB    = m0;
1388
1389     /* Optional parameters */
1390     nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1391                    &q, &m, ctypeB, &qb, &mb, check);
1392
1393     /* Nasty switch that falls thru all the way down! */
1394     if (nscan > 0)
1395     {
1396         q0 = qB = q;
1397         if (nscan > 1)
1398         {
1399             m0 = mB = m;
1400             if (nscan > 2)
1401             {
1402                 if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET)
1403                 {
1404                     auto message = gmx::formatString("Atomtype %s not found", ctypeB);
1405                     warning_error_and_exit(wi, message, FARGS);
1406                 }
1407                 qB = get_atomtype_qA(typeB, atype);
1408                 mB = get_atomtype_massA(typeB, atype);
1409                 if (nscan > 3)
1410                 {
1411                     qB = qb;
1412                     if (nscan > 4)
1413                     {
1414                         mB = mb;
1415                         if (nscan > 5)
1416                         {
1417                             warning_error(wi, "Too many parameters");
1418                         }
1419                     }
1420                 }
1421             }
1422         }
1423     }
1424
1425     push_cg(cgs, lastcg, cgnumber, nr);
1426
1427     push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype),
1428                   type, ctype, ptype, resnumberic,
1429                   resname, name, m0, q0, typeB,
1430                   typeB == type ? ctype : ctypeB, mB, qB, wi);
1431 }
1432
1433 void push_molt(t_symtab                         *symtab,
1434                std::vector<MoleculeInformation> *mol,
1435                char                             *line,
1436                warninp                          *wi)
1437 {
1438     char       type[STRLEN];
1439     int        nrexcl;
1440
1441     if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
1442     {
1443         warning_error(wi, "Expected a molecule type name and nrexcl");
1444     }
1445
1446     /* Test if this moleculetype overwrites another */
1447     const auto found = std::find_if(mol->begin(), mol->end(),
1448                                     [&type](const auto &m)
1449                                     { return strcmp(*(m.name), type) == 0; });
1450     if (found != mol->end())
1451     {
1452         auto message = gmx::formatString("moleculetype %s is redefined", type);
1453         warning_error_and_exit(wi, message, FARGS);
1454     }
1455
1456     mol->emplace_back();
1457     mol->back().initMolInfo();
1458
1459     /* Fill in the values */
1460     mol->back().name     = put_symtab(symtab, type);
1461     mol->back().nrexcl   = nrexcl;
1462     mol->back().excl_set = false;
1463 }
1464
1465 static bool default_nb_params(int ftype, gmx::ArrayRef<InteractionTypeParameters> bt, t_atoms *at,
1466                               t_param *p, int c_start, bool bB, bool bGenPairs)
1467 {
1468     int          i, j, ti, tj, ntype;
1469     bool         bFound;
1470     t_param     *pi    = nullptr;
1471     int          nr    = bt[ftype].nr;
1472     int          nral  = NRAL(ftype);
1473     int          nrfp  = interaction_function[ftype].nrfpA;
1474     int          nrfpB = interaction_function[ftype].nrfpB;
1475
1476     if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1477     {
1478         return TRUE;
1479     }
1480
1481     bFound = FALSE;
1482     if (bGenPairs)
1483     {
1484         /* First test the generated-pair position to save
1485          * time when we have 1000*1000 entries for e.g. OPLS...
1486          */
1487         ntype = static_cast<int>(std::sqrt(static_cast<double>(nr)));
1488         GMX_ASSERT(ntype * ntype == nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
1489         if (bB)
1490         {
1491             ti = at->atom[p->a[0]].typeB;
1492             tj = at->atom[p->a[1]].typeB;
1493         }
1494         else
1495         {
1496             ti = at->atom[p->a[0]].type;
1497             tj = at->atom[p->a[1]].type;
1498         }
1499         pi     = &(bt[ftype].param[ntype*ti+tj]);
1500         bFound = ((ti == pi->a[0]) && (tj == pi->a[1]));
1501     }
1502
1503     /* Search explicitly if we didnt find it */
1504     if (!bFound)
1505     {
1506         for (i = 0; ((i < nr) && !bFound); i++)
1507         {
1508             pi = &(bt[ftype].param[i]);
1509             if (bB)
1510             {
1511                 for (j = 0; ((j < nral) &&
1512                              (at->atom[p->a[j]].typeB == pi->a[j])); j++)
1513                 {
1514                     ;
1515                 }
1516             }
1517             else
1518             {
1519                 for (j = 0; ((j < nral) &&
1520                              (at->atom[p->a[j]].type == pi->a[j])); j++)
1521                 {
1522                     ;
1523                 }
1524             }
1525             bFound = (j == nral);
1526         }
1527     }
1528
1529     if (bFound)
1530     {
1531         if (bB)
1532         {
1533             if (nrfp+nrfpB > MAXFORCEPARAM)
1534             {
1535                 gmx_incons("Too many force parameters");
1536             }
1537             for (j = c_start; (j < nrfpB); j++)
1538             {
1539                 p->c[nrfp+j] = pi->c[j];
1540             }
1541         }
1542         else
1543         {
1544             for (j = c_start; (j < nrfp); j++)
1545             {
1546                 p->c[j] = pi->c[j];
1547             }
1548         }
1549     }
1550     else
1551     {
1552         for (j = c_start; (j < nrfp); j++)
1553         {
1554             p->c[j] = 0.0;
1555         }
1556     }
1557     return bFound;
1558 }
1559
1560 static bool default_cmap_params(gmx::ArrayRef<InteractionTypeParameters> bondtype,
1561                                 t_atoms *at, gpp_atomtype *atype,
1562                                 t_param *p, bool bB,
1563                                 int *cmap_type, int *nparam_def,
1564                                 warninp *wi)
1565 {
1566     int        i, nparam_found;
1567     int        ct;
1568     bool       bFound = FALSE;
1569
1570     nparam_found = 0;
1571     ct           = 0;
1572
1573     /* Match the current cmap angle against the list of cmap_types */
1574     for (i = 0; i < bondtype[F_CMAP].nct() && !bFound; i += 6)
1575     {
1576         if (bB)
1577         {
1578
1579         }
1580         else
1581         {
1582             if (
1583                 (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i])   &&
1584                 (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+1]) &&
1585                 (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+2]) &&
1586                 (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+3]) &&
1587                 (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+4]))
1588             {
1589                 /* Found cmap torsion */
1590                 bFound       = TRUE;
1591                 ct           = bondtype[F_CMAP].cmapAtomTypes[i+5];
1592                 nparam_found = 1;
1593             }
1594         }
1595     }
1596
1597     /* If we did not find a matching type for this cmap torsion */
1598     if (!bFound)
1599     {
1600         auto message = gmx::formatString("Unknown cmap torsion between atoms %d %d %d %d %d",
1601                                          p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1);
1602         warning_error_and_exit(wi, message, FARGS);
1603     }
1604
1605     *nparam_def = nparam_found;
1606     *cmap_type  = ct;
1607
1608     return bFound;
1609 }
1610
1611 /* Returns the number of exact atom type matches, i.e. non wild-card matches,
1612  * returns -1 when there are no matches at all.
1613  */
1614 static int natom_match(t_param *pi,
1615                        int type_i, int type_j, int type_k, int type_l,
1616                        const gpp_atomtype* atype)
1617 {
1618     if ((pi->ai() == -1 || get_atomtype_batype(type_i, atype) == pi->ai()) &&
1619         (pi->aj() == -1 || get_atomtype_batype(type_j, atype) == pi->aj()) &&
1620         (pi->ak() == -1 || get_atomtype_batype(type_k, atype) == pi->ak()) &&
1621         (pi->al() == -1 || get_atomtype_batype(type_l, atype) == pi->al()))
1622     {
1623         return
1624             (pi->ai() == -1 ? 0 : 1) +
1625             (pi->aj() == -1 ? 0 : 1) +
1626             (pi->ak() == -1 ? 0 : 1) +
1627             (pi->al() == -1 ? 0 : 1);
1628     }
1629     else
1630     {
1631         return -1;
1632     }
1633 }
1634
1635 static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRef<InteractionTypeParameters> bt,
1636                                              t_atoms *at, gpp_atomtype *atype,
1637                                              t_param *p, bool bB,
1638                                              t_param **param_def,
1639                                              int *nparam_def)
1640 {
1641     int          nparam_found;
1642     bool         bFound, bSame;
1643     t_param     *pi    = nullptr;
1644     t_param     *pj    = nullptr;
1645     int          nr    = bt[ftype].nr;
1646     int          nral  = NRAL(ftype);
1647     int          nrfpA = interaction_function[ftype].nrfpA;
1648     int          nrfpB = interaction_function[ftype].nrfpB;
1649
1650     if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1651     {
1652         return TRUE;
1653     }
1654
1655
1656     bFound       = FALSE;
1657     nparam_found = 0;
1658     if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1659     {
1660         int nmatch_max = -1;
1661         int i          = -1;
1662         int t;
1663
1664         /* For dihedrals we allow wildcards. We choose the first type
1665          * that has the most real matches, i.e. non-wildcard matches.
1666          */
1667         for (t = 0; ((t < nr) && nmatch_max < 4); t++)
1668         {
1669             int      nmatch;
1670             t_param *pt;
1671
1672             pt = &(bt[ftype].param[t]);
1673             if (bB)
1674             {
1675                 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);
1676             }
1677             else
1678             {
1679                 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);
1680             }
1681             if (nmatch > nmatch_max)
1682             {
1683                 nmatch_max = nmatch;
1684                 i          = t;
1685                 bFound     = TRUE;
1686             }
1687         }
1688
1689         if (bFound)
1690         {
1691             int j;
1692
1693             pi    = &(bt[ftype].param[i]);
1694             nparam_found++;
1695
1696             /* Find additional matches for this dihedral - necessary
1697              * for ftype==9.
1698              * The rule in that case is that additional matches
1699              * HAVE to be on adjacent lines!
1700              */
1701             bSame = TRUE;
1702             /* Continue from current i value */
1703             for (j = i + 2; j < nr && bSame; j += 2)
1704             {
1705                 pj    = &(bt[ftype].param[j]);
1706                 bSame = (pi->ai() == pj->ai() && pi->aj() == pj->aj() && pi->ak() == pj->ak() && pi->al() == pj->al());
1707                 if (bSame)
1708                 {
1709                     nparam_found++;
1710                 }
1711                 /* nparam_found will be increased as long as the numbers match */
1712             }
1713         }
1714     }
1715     else   /* Not a dihedral */
1716     {
1717         int i, j;
1718
1719         for (i = 0; ((i < nr) && !bFound); i++)
1720         {
1721             pi = &(bt[ftype].param[i]);
1722             if (bB)
1723             {
1724                 for (j = 0; ((j < nral) &&
1725                              (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++)
1726                 {
1727                     ;
1728                 }
1729             }
1730             else
1731             {
1732                 for (j = 0; ((j < nral) &&
1733                              (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++)
1734                 {
1735                     ;
1736                 }
1737             }
1738             bFound = (j == nral);
1739         }
1740         if (bFound)
1741         {
1742             nparam_found = 1;
1743         }
1744     }
1745
1746     *param_def  = pi;
1747     *nparam_def = nparam_found;
1748
1749     return bFound;
1750 }
1751
1752
1753
1754 void push_bond(Directive d, gmx::ArrayRef<InteractionTypeParameters> bondtype,
1755                gmx::ArrayRef<InteractionTypeParameters> bond,
1756                t_atoms *at, gpp_atomtype *atype, char *line,
1757                bool bBonded, bool bGenPairs, real fudgeQQ,
1758                bool bZero, bool *bWarn_copy_A_B,
1759                warninp *wi)
1760 {
1761     const char  *aaformat[MAXATOMLIST] = {
1762         "%d%d",
1763         "%d%d%d",
1764         "%d%d%d%d",
1765         "%d%d%d%d%d",
1766         "%d%d%d%d%d%d",
1767         "%d%d%d%d%d%d%d"
1768     };
1769     const char  *asformat[MAXATOMLIST] = {
1770         "%*s%*s",
1771         "%*s%*s%*s",
1772         "%*s%*s%*s%*s",
1773         "%*s%*s%*s%*s%*s",
1774         "%*s%*s%*s%*s%*s%*s",
1775         "%*s%*s%*s%*s%*s%*s%*s"
1776     };
1777     const char  *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1778     int          nr, i, j, nral, nral_fmt, nread, ftype;
1779     char         format[STRLEN];
1780     /* One force parameter more, so we can check if we read too many */
1781     double       cc[MAXFORCEPARAM+1];
1782     int          aa[MAXATOMLIST+1];
1783     t_param      param, *param_defA, *param_defB;
1784     bool         bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE;
1785     int          nparam_defA, nparam_defB;
1786
1787     nparam_defA = nparam_defB = 0;
1788
1789     ftype = ifunc_index(d, 1);
1790     nral  = NRAL(ftype);
1791     for (j = 0; j < MAXATOMLIST; j++)
1792     {
1793         aa[j] = NOTSET;
1794     }
1795     bDef = (NRFP(ftype) > 0);
1796
1797     if (ftype == F_SETTLE)
1798     {
1799         /* SETTLE acts on 3 atoms, but the topology format only specifies
1800          * the first atom (for historical reasons).
1801          */
1802         nral_fmt = 1;
1803     }
1804     else
1805     {
1806         nral_fmt = nral;
1807     }
1808
1809     nread = sscanf(line, aaformat[nral_fmt-1],
1810                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1811
1812     if (ftype == F_SETTLE)
1813     {
1814         aa[3] = aa[1];
1815         aa[1] = aa[0] + 1;
1816         aa[2] = aa[0] + 2;
1817     }
1818
1819     if (nread < nral_fmt)
1820     {
1821         too_few(wi);
1822         return;
1823     }
1824     else if (nread > nral_fmt)
1825     {
1826         /* this is a hack to allow for virtual sites with swapped parity */
1827         bSwapParity = (aa[nral] < 0);
1828         if (bSwapParity)
1829         {
1830             aa[nral] = -aa[nral];
1831         }
1832         ftype = ifunc_index(d, aa[nral]);
1833         if (bSwapParity)
1834         {
1835             switch (ftype)
1836             {
1837                 case F_VSITE3FAD:
1838                 case F_VSITE3OUT:
1839                     break;
1840                 default:
1841                     auto message = gmx::formatString("Negative function types only allowed for %s and %s",
1842                                                      interaction_function[F_VSITE3FAD].longname,
1843                                                      interaction_function[F_VSITE3OUT].longname);
1844                     warning_error_and_exit(wi, message, FARGS);
1845             }
1846         }
1847     }
1848
1849
1850     /* Check for double atoms and atoms out of bounds */
1851     for (i = 0; (i < nral); i++)
1852     {
1853         if (aa[i] < 1 || aa[i] > at->nr)
1854         {
1855             auto message = gmx::formatString
1856                     ("Atom index (%d) in %s out of bounds (1-%d).\n"
1857                     "This probably means that you have inserted topology section \"%s\"\n"
1858                     "in a part belonging to a different molecule than you intended to.\n"
1859                     "In that case move the \"%s\" section to the right molecule.",
1860                     aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
1861             warning_error_and_exit(wi, message, FARGS);
1862         }
1863         for (j = i+1; (j < nral); j++)
1864         {
1865             GMX_ASSERT(j < MAXATOMLIST + 1, "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1866             if (aa[i] == aa[j])
1867             {
1868                 auto message = gmx::formatString("Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1869                 if (ftype == F_ANGRES)
1870                 {
1871                     /* Since the angle restraints uses 2 pairs of atoms to
1872                      * defines an angle between vectors, it can be useful
1873                      * to use one atom twice, so we only issue a note here.
1874                      */
1875                     warning_note(wi, message);
1876                 }
1877                 else
1878                 {
1879                     warning_error(wi, message);
1880                 }
1881             }
1882         }
1883     }
1884
1885     /* default force parameters  */
1886     for (j = 0; (j < MAXATOMLIST); j++)
1887     {
1888         param.a[j] = aa[j]-1;
1889     }
1890     for (j = 0; (j < MAXFORCEPARAM); j++)
1891     {
1892         param.c[j] = 0.0;
1893     }
1894
1895     /* Get force params for normal and free energy perturbation
1896      * studies, as determined by types!
1897      */
1898
1899     if (bBonded)
1900     {
1901         bFoundA = defaultInteractionTypeParameters(ftype, bondtype, at, atype, &param, FALSE, &param_defA, &nparam_defA);
1902         if (bFoundA)
1903         {
1904             /* Copy the A-state and B-state default parameters. */
1905             GMX_ASSERT(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM, "Bonded interactions may have at most 12 parameters");
1906             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
1907             {
1908                 param.c[j] = param_defA->c[j];
1909             }
1910         }
1911         bFoundB = defaultInteractionTypeParameters(ftype, bondtype, at, atype, &param, TRUE, &param_defB, &nparam_defB);
1912         if (bFoundB)
1913         {
1914             /* Copy only the B-state default parameters */
1915             for (j = NRFPA(ftype); (j < NRFP(ftype)); j++)
1916             {
1917                 param.c[j] = param_defB->c[j];
1918             }
1919         }
1920     }
1921     else if (ftype == F_LJ14)
1922     {
1923         bFoundA = default_nb_params(ftype, bondtype, at, &param, 0, FALSE, bGenPairs);
1924         bFoundB = default_nb_params(ftype, bondtype, at, &param, 0, TRUE, bGenPairs);
1925     }
1926     else if (ftype == F_LJC14_Q)
1927     {
1928         param.c[0] = fudgeQQ;
1929         /* Fill in the A-state charges as default parameters */
1930         param.c[1] = at->atom[param.a[0]].q;
1931         param.c[2] = at->atom[param.a[1]].q;
1932         /* The default LJ parameters are the standard 1-4 parameters */
1933         bFoundA = default_nb_params(F_LJ14, bondtype, at, &param, 3, FALSE, bGenPairs);
1934         bFoundB = TRUE;
1935     }
1936     else if (ftype == F_LJC_PAIRS_NB)
1937     {
1938         /* Defaults are not supported here */
1939         bFoundA = FALSE;
1940         bFoundB = TRUE;
1941     }
1942     else
1943     {
1944         gmx_incons("Unknown function type in push_bond");
1945     }
1946
1947     if (nread > nral_fmt)
1948     {
1949         /* Manually specified parameters - in this case we discard multiple torsion info! */
1950
1951         strcpy(format, asformat[nral_fmt-1]);
1952         strcat(format, ccformat);
1953
1954         nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5],
1955                        &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
1956
1957         if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1958         {
1959             /* We only have to issue a warning if these atoms are perturbed! */
1960             bPert = FALSE;
1961             for (j = 0; (j < nral); j++)
1962             {
1963                 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1964             }
1965
1966             if (bPert && *bWarn_copy_A_B)
1967             {
1968                 auto message = gmx::formatString("Some parameters for bonded interaction involving "
1969                                                  "perturbed atoms are specified explicitly in "
1970                                                  "state A, but not B - copying A to B");
1971                 warning(wi, message);
1972                 *bWarn_copy_A_B = FALSE;
1973             }
1974
1975             /* If only the A parameters were specified, copy them to the B state */
1976             /* The B-state parameters correspond to the first nrfpB
1977              * A-state parameters.
1978              */
1979             for (j = 0; (j < NRFPB(ftype)); j++)
1980             {
1981                 cc[nread++] = cc[j];
1982             }
1983         }
1984
1985         /* If nread was 0 or EOF, no parameters were read => use defaults.
1986          * If nread was nrfpA we copied above so nread=nrfp.
1987          * If nread was nrfp we are cool.
1988          * For F_LJC14_Q we allow supplying fudgeQQ only.
1989          * Anything else is an error!
1990          */
1991         if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1992             !(ftype == F_LJC14_Q && nread == 1))
1993         {
1994             auto message = gmx::formatString
1995                     ("Incorrect number of parameters - found %d, expected %d "
1996                     "or %d for %s (after the function type).",
1997                     nread, NRFPA(ftype), NRFP(ftype), interaction_function[ftype].longname);
1998             warning_error_and_exit(wi, message, FARGS);
1999         }
2000
2001         for (j = 0; (j < nread); j++)
2002         {
2003             param.c[j] = cc[j];
2004         }
2005
2006         /* Check whether we have to use the defaults */
2007         if (nread == NRFP(ftype))
2008         {
2009             bDef = FALSE;
2010         }
2011     }
2012     else
2013     {
2014         nread = 0;
2015     }
2016     /* nread now holds the number of force parameters read! */
2017
2018     if (bDef)
2019     {
2020         /* Use defaults */
2021         /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
2022         if (ftype == F_PDIHS)
2023         {
2024             if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB)))
2025             {
2026                 auto message =
2027                     gmx::formatString("Cannot automatically perturb a torsion with multiple terms to different form.\n"
2028                                       "Please specify perturbed parameters manually for this torsion in your topology!");
2029                 warning_error(wi, message);
2030             }
2031         }
2032
2033         if (nread > 0 && nread < NRFPA(ftype))
2034         {
2035             /* Issue an error, do not use defaults */
2036             auto message = gmx::formatString("Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
2037             warning_error(wi, message);
2038         }
2039
2040         if (nread == 0 || nread == EOF)
2041         {
2042             if (!bFoundA)
2043             {
2044                 if (interaction_function[ftype].flags & IF_VSITE)
2045                 {
2046                     /* set them to NOTSET, will be calculated later */
2047                     for (j = 0; (j < MAXFORCEPARAM); j++)
2048                     {
2049                         param.c[j] = NOTSET;
2050                     }
2051
2052                     if (bSwapParity)
2053                     {
2054                         param.c1() = -1; /* flag to swap parity of vsite construction */
2055                     }
2056                 }
2057                 else
2058                 {
2059                     if (bZero)
2060                     {
2061                         fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
2062                                 interaction_function[ftype].longname);
2063                     }
2064                     else
2065                     {
2066                         auto message = gmx::formatString("No default %s types", interaction_function[ftype].longname);
2067                         warning_error(wi, message);
2068                     }
2069                 }
2070             }
2071             else
2072             {
2073                 if (bSwapParity)
2074                 {
2075                     switch (ftype)
2076                     {
2077                         case F_VSITE3FAD:
2078                             param.c0() = 360 - param.c0();
2079                             break;
2080                         case F_VSITE3OUT:
2081                             param.c2() = -param.c2();
2082                             break;
2083                     }
2084                 }
2085             }
2086             if (!bFoundB)
2087             {
2088                 /* We only have to issue a warning if these atoms are perturbed! */
2089                 bPert = FALSE;
2090                 for (j = 0; (j < nral); j++)
2091                 {
2092                     bPert = bPert || PERTURBED(at->atom[param.a[j]]);
2093                 }
2094
2095                 if (bPert)
2096                 {
2097                     auto message = gmx::formatString("No default %s types for perturbed atoms, "
2098                                                      "using normal values", interaction_function[ftype].longname);
2099                     warning(wi, message);
2100                 }
2101             }
2102         }
2103     }
2104
2105     if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
2106         && param.c[5] != param.c[2])
2107     {
2108         auto message = gmx::formatString("%s multiplicity can not be perturbed %f!=%f",
2109                                          interaction_function[ftype].longname,
2110                                          param.c[2], param.c[5]);
2111         warning_error_and_exit(wi, message, FARGS);
2112     }
2113
2114     if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
2115     {
2116         auto message = gmx::formatString("%s table number can not be perturbed %d!=%d",
2117                                          interaction_function[ftype].longname,
2118                                          gmx::roundToInt(param.c[0]), gmx::roundToInt(param.c[2]));
2119         warning_error_and_exit(wi, message, FARGS);
2120     }
2121
2122     /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2123     if (ftype == F_RBDIHS)
2124     {
2125         nr = 0;
2126         for (i = 0; i < NRFP(ftype); i++)
2127         {
2128             if (param.c[i] != 0)
2129             {
2130                 nr++;
2131             }
2132         }
2133         if (nr == 0)
2134         {
2135             return;
2136         }
2137     }
2138
2139     /* Put the values in the appropriate arrays */
2140     add_param_to_list (&bond[ftype], &param);
2141
2142     /* Push additional torsions from FF for ftype==9 if we have them.
2143      * We have already checked that the A/B states do not differ in this case,
2144      * so we do not have to double-check that again, or the vsite stuff.
2145      * In addition, those torsions cannot be automatically perturbed.
2146      */
2147     if (bDef && ftype == F_PDIHS)
2148     {
2149         for (i = 1; i < nparam_defA; i++)
2150         {
2151             /* Advance pointer! */
2152             param_defA += 2;
2153             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
2154             {
2155                 param.c[j] = param_defA->c[j];
2156             }
2157             /* And push the next term for this torsion */
2158             add_param_to_list (&bond[ftype], &param);
2159         }
2160     }
2161 }
2162
2163 void push_cmap(Directive d, gmx::ArrayRef<InteractionTypeParameters> bondtype,
2164                gmx::ArrayRef<InteractionTypeParameters> bond,
2165                t_atoms *at, gpp_atomtype *atype, char *line,
2166                warninp *wi)
2167 {
2168     const char *aaformat[MAXATOMLIST+1] =
2169     {
2170         "%d",
2171         "%d%d",
2172         "%d%d%d",
2173         "%d%d%d%d",
2174         "%d%d%d%d%d",
2175         "%d%d%d%d%d%d",
2176         "%d%d%d%d%d%d%d"
2177     };
2178
2179     int      i, j, ftype, nral, nread, ncmap_params;
2180     int      cmap_type;
2181     int      aa[MAXATOMLIST+1];
2182     bool     bFound;
2183     t_param  param;
2184
2185     ftype        = ifunc_index(d, 1);
2186     nral         = NRAL(ftype);
2187     ncmap_params = 0;
2188
2189     nread = sscanf(line, aaformat[nral-1],
2190                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2191
2192     if (nread < nral)
2193     {
2194         too_few(wi);
2195         return;
2196     }
2197     else if (nread == nral)
2198     {
2199         ftype = ifunc_index(d, 1);
2200     }
2201
2202     /* Check for double atoms and atoms out of bounds */
2203     for (i = 0; i < nral; i++)
2204     {
2205         if (aa[i] < 1 || aa[i] > at->nr)
2206         {
2207             auto message = gmx::formatString
2208                     ("Atom index (%d) in %s out of bounds (1-%d).\n"
2209                     "This probably means that you have inserted topology section \"%s\"\n"
2210                     "in a part belonging to a different molecule than you intended to.\n"
2211                     "In that case move the \"%s\" section to the right molecule.",
2212                     aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2213             warning_error_and_exit(wi, message, FARGS);
2214         }
2215
2216         for (j = i+1; (j < nral); j++)
2217         {
2218             if (aa[i] == aa[j])
2219             {
2220                 auto message = gmx::formatString("Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2221                 warning_error(wi, message);
2222             }
2223         }
2224     }
2225
2226     /* default force parameters  */
2227     for (j = 0; (j < MAXATOMLIST); j++)
2228     {
2229         param.a[j] = aa[j]-1;
2230     }
2231     for (j = 0; (j < MAXFORCEPARAM); j++)
2232     {
2233         param.c[j] = 0.0;
2234     }
2235
2236     /* Get the cmap type for this cmap angle */
2237     bFound = default_cmap_params(bondtype, at, atype, &param, FALSE, &cmap_type, &ncmap_params, wi);
2238
2239     /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2240     if (bFound && ncmap_params == 1)
2241     {
2242         /* Put the values in the appropriate arrays */
2243         param.c[0] = cmap_type;
2244         add_param_to_list(&bond[ftype], &param);
2245     }
2246     else
2247     {
2248         /* This is essentially the same check as in default_cmap_params() done one more time */
2249         auto message = gmx::formatString("Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2250                                          param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2251         warning_error_and_exit(wi, message, FARGS);
2252     }
2253 }
2254
2255
2256
2257 void push_vsitesn(Directive d, gmx::ArrayRef<InteractionTypeParameters> bond,
2258                   t_atoms *at, char *line,
2259                   warninp *wi)
2260 {
2261     char   *ptr;
2262     int     type, ftype, j, n, ret, nj, a;
2263     int    *atc    = nullptr;
2264     double *weight = nullptr, weight_tot;
2265     t_param param;
2266
2267     /* default force parameters  */
2268     for (j = 0; (j < MAXATOMLIST); j++)
2269     {
2270         param.a[j] = NOTSET;
2271     }
2272     for (j = 0; (j < MAXFORCEPARAM); j++)
2273     {
2274         param.c[j] = 0.0;
2275     }
2276
2277     ptr  = line;
2278     ret  = sscanf(ptr, "%d%n", &a, &n);
2279     ptr += n;
2280     if (ret == 0)
2281     {
2282         auto message = gmx::formatString("Expected an atom index in section \"%s\"",
2283                                          dir2str(d));
2284         warning_error_and_exit(wi, message, FARGS);
2285     }
2286
2287     param.a[0] = a - 1;
2288
2289     sscanf(ptr, "%d%n", &type, &n);
2290     ptr  += n;
2291     ftype = ifunc_index(d, type);
2292
2293     weight_tot = 0;
2294     nj         = 0;
2295     do
2296     {
2297         ret  = sscanf(ptr, "%d%n", &a, &n);
2298         ptr += n;
2299         if (ret > 0)
2300         {
2301             if (nj % 20 == 0)
2302             {
2303                 srenew(atc, nj+20);
2304                 srenew(weight, nj+20);
2305             }
2306             atc[nj] = a - 1;
2307             switch (type)
2308             {
2309                 case 1:
2310                     weight[nj] = 1;
2311                     break;
2312                 case 2:
2313                     /* Here we use the A-state mass as a parameter.
2314                      * Note that the B-state mass has no influence.
2315                      */
2316                     weight[nj] = at->atom[atc[nj]].m;
2317                     break;
2318                 case 3:
2319                     weight[nj] = -1;
2320                     ret        = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2321                     ptr       += n;
2322                     if (weight[nj] < 0)
2323                     {
2324                         auto message = gmx::formatString
2325                                 ("No weight or negative weight found for vsiten "
2326                                 "constructing atom %d (atom index %d)",
2327                                 nj+1, atc[nj]+1);
2328                         warning_error_and_exit(wi, message, FARGS);
2329                     }
2330                     break;
2331                 default:
2332                     auto message = gmx::formatString("Unknown vsiten type %d", type);
2333                     warning_error_and_exit(wi, message, FARGS);
2334             }
2335             weight_tot += weight[nj];
2336             nj++;
2337         }
2338     }
2339     while (ret > 0);
2340
2341     if (nj == 0)
2342     {
2343         auto message = gmx::formatString("Expected more than one atom index in section \"%s\"",
2344                                          dir2str(d));
2345         warning_error_and_exit(wi, message, FARGS);
2346     }
2347
2348     if (weight_tot == 0)
2349     {
2350         warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
2351     }
2352
2353     for (j = 0; j < nj; j++)
2354     {
2355         param.a[1] = atc[j];
2356         param.c[0] = nj;
2357         param.c[1] = weight[j]/weight_tot;
2358         /* Put the values in the appropriate arrays */
2359         add_param_to_list (&bond[ftype], &param);
2360     }
2361
2362     sfree(atc);
2363     sfree(weight);
2364 }
2365
2366 void push_mol(gmx::ArrayRef<MoleculeInformation> mols, char *pline, int *whichmol,
2367               int *nrcopies,
2368               warninp *wi)
2369 {
2370     char type[STRLEN];
2371
2372     if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2373     {
2374         too_few(wi);
2375         return;
2376     }
2377
2378     /* Search moleculename.
2379      * Here we originally only did case insensitive matching. But because
2380      * some PDB files can have many chains and use case to generate more
2381      * chain-identifiers, which in turn end up in our moleculetype name,
2382      * we added support for case-sensitivity.
2383      */
2384     int nrcs    = 0;
2385     int nrci    = 0;
2386     int matchci = -1;
2387     int matchcs = -1;
2388     int i       = 0;
2389     for (const auto &mol : mols)
2390     {
2391         if (strcmp(type, *(mol.name)) == 0)
2392         {
2393             nrcs++;
2394             matchcs = i;
2395         }
2396         if (gmx_strcasecmp(type, *(mol.name)) == 0)
2397         {
2398             nrci++;
2399             matchci = i;
2400         }
2401         i++;
2402     }
2403
2404     if (nrcs == 1)
2405     {
2406         // select the case sensitive match
2407         *whichmol = matchcs;
2408     }
2409     else
2410     {
2411         // avoid matching case-insensitive when we have multiple matches
2412         if (nrci > 1)
2413         {
2414             auto message = gmx::formatString
2415                     ("For moleculetype '%s' in [ system ] %d case insensitive "
2416                     "matches, but %d case sensitive matches were found. Check "
2417                     "the case of the characters in the moleculetypes.",
2418                     type, nrci, nrcs);
2419             warning_error_and_exit(wi, message, FARGS);
2420         }
2421         if (nrci == 1)
2422         {
2423             // select the unique case insensitive match
2424             *whichmol = matchci;
2425         }
2426         else
2427         {
2428             auto message = gmx::formatString("No such moleculetype %s", type);
2429             warning_error_and_exit(wi, message, FARGS);
2430         }
2431     }
2432 }
2433
2434 void push_excl(char *line, gmx::ArrayRef<gmx::ExclusionBlock> b2, warninp *wi)
2435 {
2436     int  i, j;
2437     int  n;
2438     char base[STRLEN], format[STRLEN];
2439
2440     if (sscanf(line, "%d", &i) == 0)
2441     {
2442         return;
2443     }
2444
2445     if ((1 <= i) && (i <= b2.ssize()))
2446     {
2447         i--;
2448     }
2449     else
2450     {
2451         return;
2452     }
2453     strcpy(base, "%*d");
2454     do
2455     {
2456         strcpy(format, base);
2457         strcat(format, "%d");
2458         n = sscanf(line, format, &j);
2459         if (n == 1)
2460         {
2461             if ((1 <= j) && (j <= b2.ssize()))
2462             {
2463                 j--;
2464                 b2[i].atomNumber.push_back(j);
2465                 /* also add the reverse exclusion! */
2466                 b2[j].atomNumber.push_back(i);
2467                 strcat(base, "%*d");
2468             }
2469             else
2470             {
2471                 auto message = gmx::formatString("Invalid Atomnr j: %d, b2->nr: %zu\n", j, b2.size());
2472                 warning_error_and_exit(wi, message, FARGS);
2473             }
2474         }
2475     }
2476     while (n == 1);
2477 }
2478
2479 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype *at,
2480                            t_nbparam ***nbparam, t_nbparam ***pair)
2481 {
2482     t_atom  atom;
2483     t_param param;
2484     int     i, nr;
2485
2486     /* Define an atom type with all parameters set to zero (no interactions) */
2487     atom.q = 0.0;
2488     atom.m = 0.0;
2489     /* Type for decoupled atoms could be anything,
2490      * this should be changed automatically later when required.
2491      */
2492     atom.ptype = eptAtom;
2493     for (i = 0; (i < MAXFORCEPARAM); i++)
2494     {
2495         param.c[i] = 0.0;
2496     }
2497
2498     nr = add_atomtype(at, symtab, &atom, "decoupled", &param, -1, 0);
2499
2500     /* Add space in the non-bonded parameters matrix */
2501     realloc_nb_params(at, nbparam, pair);
2502
2503     return nr;
2504 }
2505
2506 static void convert_pairs_to_pairsQ(gmx::ArrayRef<InteractionTypeParameters> plist,
2507                                     real fudgeQQ, t_atoms *atoms)
2508 {
2509     t_param *paramp1, *paramp2, *paramnew;
2510     int      i, j, p1nr, p2nr, p2newnr;
2511
2512     /* Add the pair list to the pairQ list */
2513     p1nr    = plist[F_LJ14].nr;
2514     p2nr    = plist[F_LJC14_Q].nr;
2515     p2newnr = p1nr + p2nr;
2516     snew(paramnew, p2newnr);
2517
2518     paramp1             = plist[F_LJ14].param;
2519     paramp2             = plist[F_LJC14_Q].param;
2520
2521     /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2522        it may be possible to just ADD the converted F_LJ14 array
2523        to the old F_LJC14_Q array, but since we have to create
2524        a new sized memory structure, better just to deep copy it all.
2525      */
2526
2527     for (i = 0; i < p2nr; i++)
2528     {
2529         /* Copy over parameters */
2530         for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */
2531         {
2532             paramnew[i].c[j] = paramp2[i].c[j];
2533         }
2534
2535         /* copy over atoms */
2536         for (j = 0; j < 2; j++)
2537         {
2538             paramnew[i].a[j] = paramp2[i].a[j];
2539         }
2540     }
2541
2542     for (i = p2nr; i < p2newnr; i++)
2543     {
2544         j             = i-p2nr;
2545
2546         /* Copy over parameters */
2547         paramnew[i].c[0] = fudgeQQ;
2548         paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q;
2549         paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q;
2550         paramnew[i].c[3] = paramp1[j].c[0];
2551         paramnew[i].c[4] = paramp1[j].c[1];
2552
2553         /* copy over atoms */
2554         paramnew[i].a[0] = paramp1[j].a[0];
2555         paramnew[i].a[1] = paramp1[j].a[1];
2556     }
2557
2558     /* free the old pairlists */
2559     sfree(plist[F_LJC14_Q].param);
2560     sfree(plist[F_LJ14].param);
2561
2562     /* now assign the new data to the F_LJC14_Q structure */
2563     plist[F_LJC14_Q].param   = paramnew;
2564     plist[F_LJC14_Q].nr      = p2newnr;
2565
2566     /* Empty the LJ14 pairlist */
2567     plist[F_LJ14].nr    = 0;
2568     plist[F_LJ14].param = nullptr;
2569 }
2570
2571 static void generate_LJCpairsNB(MoleculeInformation *mol, int nb_funct, InteractionTypeParameters *nbp, warninp *wi)
2572 {
2573     int       n, ntype, i, j, k;
2574     t_atom   *atom;
2575     t_blocka *excl;
2576     bool      bExcl;
2577     t_param   param;
2578
2579     n    = mol->atoms.nr;
2580     atom = mol->atoms.atom;
2581
2582     ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->nr)));
2583     GMX_ASSERT(ntype * ntype == nbp->nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
2584
2585     for (i = 0; i < MAXATOMLIST; i++)
2586     {
2587         param.a[i] = NOTSET;
2588     }
2589     for (i = 0; i < MAXFORCEPARAM; i++)
2590     {
2591         param.c[i] = NOTSET;
2592     }
2593
2594     /* Add a pair interaction for all non-excluded atom pairs */
2595     excl = &mol->excls;
2596     for (i = 0; i < n; i++)
2597     {
2598         for (j = i+1; j < n; j++)
2599         {
2600             bExcl = FALSE;
2601             for (k = excl->index[i]; k < excl->index[i+1]; k++)
2602             {
2603                 if (excl->a[k] == j)
2604                 {
2605                     bExcl = TRUE;
2606                 }
2607             }
2608             if (!bExcl)
2609             {
2610                 if (nb_funct != F_LJ)
2611                 {
2612                     auto message = gmx::formatString
2613                             ("Can only generate non-bonded pair interactions "
2614                             "for Van der Waals type Lennard-Jones");
2615                     warning_error_and_exit(wi, message, FARGS);
2616                 }
2617                 param.a[0] = i;
2618                 param.a[1] = j;
2619                 param.c[0] = atom[i].q;
2620                 param.c[1] = atom[j].q;
2621                 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2622                 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2623                 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], &param);
2624             }
2625         }
2626     }
2627 }
2628
2629 static void set_excl_all(t_blocka *excl)
2630 {
2631     int nat, i, j, k;
2632
2633     /* Get rid of the current exclusions and exclude all atom pairs */
2634     nat       = excl->nr;
2635     excl->nra = nat*nat;
2636     srenew(excl->a, excl->nra);
2637     k = 0;
2638     for (i = 0; i < nat; i++)
2639     {
2640         excl->index[i] = k;
2641         for (j = 0; j < nat; j++)
2642         {
2643             excl->a[k++] = j;
2644         }
2645     }
2646     excl->index[nat] = k;
2647 }
2648
2649 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2650                            int couple_lam0, int couple_lam1,
2651                            const char *mol_name, warninp *wi)
2652 {
2653     int  i;
2654
2655     for (i = 0; i < atoms->nr; i++)
2656     {
2657         t_atom *atom;
2658
2659         atom = &atoms->atom[i];
2660
2661         if (atom->qB != atom->q || atom->typeB != atom->type)
2662         {
2663             auto message = gmx::formatString
2664                     ("Atom %d in molecule type '%s' has different A and B state "
2665                     "charges and/or atom types set in the topology file as well "
2666                     "as through the mdp option '%s'. You can not use both "
2667                     "these methods simultaneously.",
2668                     i + 1, mol_name, "couple-moltype");
2669             warning_error_and_exit(wi, message, FARGS);
2670         }
2671
2672         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2673         {
2674             atom->q     = 0.0;
2675         }
2676         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2677         {
2678             atom->type  = atomtype_decouple;
2679         }
2680         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2681         {
2682             atom->qB    = 0.0;
2683         }
2684         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2685         {
2686             atom->typeB = atomtype_decouple;
2687         }
2688     }
2689 }
2690
2691 void convert_moltype_couple(MoleculeInformation *mol, int atomtype_decouple, real fudgeQQ,
2692                             int couple_lam0, int couple_lam1,
2693                             bool bCoupleIntra, int nb_funct, InteractionTypeParameters *nbp,
2694                             warninp *wi)
2695 {
2696     convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2697     if (!bCoupleIntra)
2698     {
2699         generate_LJCpairsNB(mol, nb_funct, nbp, wi);
2700         set_excl_all(&mol->excls);
2701     }
2702     decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1,
2703                    *mol->name, wi);
2704 }