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