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