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