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     const gmx::ArrayRef<const int> bondTypeAtomTypes = bondType.atoms();
1731     GMX_RELEASE_ASSERT(bondTypeAtomTypes.size() == 4, "Dihedral types have 4 atom types");
1732     int numExactMatches = 0;
1733     if (std::equal(bondTypeAtomTypes.begin(),
1734                    bondTypeAtomTypes.end(),
1735                    atomTypes.begin(),
1736                    [&numExactMatches](int bondTypeAtomType, int atomType) {
1737                        if (bondTypeAtomType == atomType)
1738                        {
1739                            // Found an exact atom type match
1740                            ++numExactMatches;
1741                            return true;
1742                        }
1743                        else if (bondTypeAtomType == -1)
1744                        {
1745                            // Found a wildcard atom type match
1746                            return true;
1747                        }
1748                        // Atom types do not match
1749                        return false;
1750                    }))
1751     {
1752         return numExactMatches;
1753     }
1754     return -1;
1755 }
1756
1757 static std::vector<InteractionOfType>::iterator
1758 defaultInteractionsOfType(int                               ftype,
1759                           gmx::ArrayRef<InteractionsOfType> bondType,
1760                           const gmx::ArrayRef<const int>    atomTypes,
1761                           int*                              nparam_def)
1762 {
1763     int nparam_found = 0;
1764
1765     if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1766     {
1767         int nmatch_max = -1;
1768
1769         /* For dihedrals we allow wildcards. We choose the first type
1770          * that has the most exact matches, i.e. non-wildcard matches.
1771          */
1772         auto prevPos = bondType[ftype].interactionTypes.end();
1773         auto pos     = bondType[ftype].interactionTypes.begin();
1774         while (pos != bondType[ftype].interactionTypes.end() && nmatch_max < 4)
1775         {
1776             pos = std::find_if(bondType[ftype].interactionTypes.begin(),
1777                                bondType[ftype].interactionTypes.end(),
1778                                [&atomTypes, &nmatch_max](const auto& bondType) {
1779                                    return (findNumberOfDihedralAtomMatches(bondType, atomTypes) > nmatch_max);
1780                                });
1781             if (pos != bondType[ftype].interactionTypes.end())
1782             {
1783                 prevPos    = pos;
1784                 nmatch_max = findNumberOfDihedralAtomMatches(*pos, atomTypes);
1785             }
1786         }
1787
1788         if (prevPos != bondType[ftype].interactionTypes.end())
1789         {
1790             nparam_found++;
1791
1792             /* Find additional matches for this dihedral - necessary
1793              * for ftype==9.
1794              * The rule in that case is that additional matches
1795              * HAVE to be on adjacent lines!
1796              */
1797             bool bSame = true;
1798             // Advance iterator (like std::advance) without incrementing past end (UB)
1799             const auto safeAdvance = [](auto& it, auto n, auto end) {
1800                 it = end - it > n ? it + n : end;
1801             };
1802             /* Continue from current iterator position */
1803             auto       nextPos = prevPos;
1804             const auto endIter = bondType[ftype].interactionTypes.end();
1805             safeAdvance(nextPos, 2, endIter);
1806             for (; nextPos < endIter && bSame; safeAdvance(nextPos, 2, endIter))
1807             {
1808                 bSame = (prevPos->ai() == nextPos->ai() && prevPos->aj() == nextPos->aj()
1809                          && prevPos->ak() == nextPos->ak() && prevPos->al() == nextPos->al());
1810                 if (bSame)
1811                 {
1812                     nparam_found++;
1813                 }
1814                 /* nparam_found will be increased as long as the numbers match */
1815             }
1816         }
1817         *nparam_def = nparam_found;
1818         return prevPos;
1819     }
1820     else /* Not a dihedral */
1821     {
1822         auto found = std::find_if(
1823                 bondType[ftype].interactionTypes.begin(),
1824                 bondType[ftype].interactionTypes.end(),
1825                 [&atomTypes](const auto& param) {
1826                     return std::equal(param.atoms().begin(), param.atoms().end(), atomTypes.begin());
1827                 });
1828         if (found != bondType[ftype].interactionTypes.end())
1829         {
1830             nparam_found = 1;
1831         }
1832         *nparam_def = nparam_found;
1833         return found;
1834     }
1835 }
1836
1837
1838 void push_bond(Directive                         d,
1839                gmx::ArrayRef<InteractionsOfType> bondtype,
1840                gmx::ArrayRef<InteractionsOfType> bond,
1841                t_atoms*                          at,
1842                PreprocessingAtomTypes*           atypes,
1843                char*                             line,
1844                bool                              bBonded,
1845                bool                              bGenPairs,
1846                real                              fudgeQQ,
1847                bool                              bZero,
1848                bool*                             bWarn_copy_A_B,
1849                warninp*                          wi)
1850 {
1851     const char* aaformat[MAXATOMLIST] = { "%d%d",       "%d%d%d",       "%d%d%d%d",
1852                                           "%d%d%d%d%d", "%d%d%d%d%d%d", "%d%d%d%d%d%d%d" };
1853     const char* asformat[MAXATOMLIST] = {
1854         "%*s%*s",          "%*s%*s%*s",          "%*s%*s%*s%*s",
1855         "%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s%*s"
1856     };
1857     const char* ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1858     int         nral, nral_fmt, nread, ftype;
1859     char        format[STRLEN];
1860     /* One force parameter more, so we can check if we read too many */
1861     double                           cc[MAXFORCEPARAM + 1];
1862     std::array<int, MAXATOMLIST + 1> aa;
1863     bool                             bFoundA = FALSE, bFoundB = FALSE, bDef, bSwapParity = FALSE;
1864     int                              nparam_defA, nparam_defB;
1865
1866     nparam_defA = nparam_defB = 0;
1867
1868     ftype = ifunc_index(d, 1);
1869     nral  = NRAL(ftype);
1870     for (int j = 0; j < nral; j++)
1871     {
1872         aa[j] = NOTSET;
1873     }
1874     bDef = (NRFP(ftype) > 0);
1875
1876     if (ftype == F_SETTLE)
1877     {
1878         /* SETTLE acts on 3 atoms, but the topology format only specifies
1879          * the first atom (for historical reasons).
1880          */
1881         nral_fmt = 1;
1882     }
1883     else
1884     {
1885         nral_fmt = nral;
1886     }
1887
1888     nread = sscanf(line, aaformat[nral_fmt - 1], &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1889
1890     if (ftype == F_SETTLE)
1891     {
1892         aa[3] = aa[1];
1893         aa[1] = aa[0] + 1;
1894         aa[2] = aa[0] + 2;
1895     }
1896
1897     if (nread < nral_fmt)
1898     {
1899         too_few(wi);
1900         return;
1901     }
1902     else if (nread > nral_fmt)
1903     {
1904         /* this is a hack to allow for virtual sites with swapped parity */
1905         bSwapParity = (aa[nral] < 0);
1906         if (bSwapParity)
1907         {
1908             aa[nral] = -aa[nral];
1909         }
1910         ftype = ifunc_index(d, aa[nral]);
1911         if (bSwapParity)
1912         {
1913             switch (ftype)
1914             {
1915                 case F_VSITE3FAD:
1916                 case F_VSITE3OUT: break;
1917                 default:
1918                     auto message =
1919                             gmx::formatString("Negative function types only allowed for %s and %s",
1920                                               interaction_function[F_VSITE3FAD].longname,
1921                                               interaction_function[F_VSITE3OUT].longname);
1922                     warning_error_and_exit(wi, message, FARGS);
1923             }
1924         }
1925     }
1926
1927
1928     /* Check for double atoms and atoms out of bounds, then convert to 0-based indexing */
1929     for (int i = 0; (i < nral); i++)
1930     {
1931         if (aa[i] < 1 || aa[i] > at->nr)
1932         {
1933             auto message = gmx::formatString(
1934                     "Atom index (%d) in %s out of bounds (1-%d).\n"
1935                     "This probably means that you have inserted topology section \"%s\"\n"
1936                     "in a part belonging to a different molecule than you intended to.\n"
1937                     "In that case move the \"%s\" section to the right molecule.",
1938                     aa[i],
1939                     enumValueToString(d),
1940                     at->nr,
1941                     enumValueToString(d),
1942                     enumValueToString(d));
1943             warning_error_and_exit(wi, message, FARGS);
1944         }
1945         for (int j = i + 1; (j < nral); j++)
1946         {
1947             GMX_ASSERT(j < MAXATOMLIST + 1,
1948                        "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1949             if (aa[i] == aa[j])
1950             {
1951                 auto message = gmx::formatString(
1952                         "Duplicate atom index (%d) in %s", aa[i], enumValueToString(d));
1953                 if (ftype == F_ANGRES)
1954                 {
1955                     /* Since the angle restraints uses 2 pairs of atoms to
1956                      * defines an angle between vectors, it can be useful
1957                      * to use one atom twice, so we only issue a note here.
1958                      */
1959                     warning_note(wi, message);
1960                 }
1961                 else
1962                 {
1963                     warning_error(wi, message);
1964                 }
1965             }
1966         }
1967
1968         // Convert to 0-based indexing
1969         --aa[i];
1970     }
1971
1972     // These are the atom indices for this interaction
1973     gmx::ArrayRef<int> atomIndices(aa.begin(), aa.begin() + nral);
1974
1975     // Look up the A-state atom types for this interaction
1976     std::vector<int> atomTypes(atomIndices.size());
1977     std::transform(atomIndices.begin(), atomIndices.end(), atomTypes.begin(), [at, atypes](const int atomIndex) {
1978         return atypes->bondAtomTypeFromAtomType(at->atom[atomIndex].type).value();
1979     });
1980     // Look up the B-state atom types for this interaction
1981     std::vector<int> atomTypesB(atomIndices.size());
1982     std::transform(atomIndices.begin(), atomIndices.end(), atomTypesB.begin(), [at, atypes](const int atomIndex) {
1983         return atypes->bondAtomTypeFromAtomType(at->atom[atomIndex].typeB).value();
1984     });
1985
1986     /* default force parameters  */
1987     /* need to have an empty but initialized param array for some reason */
1988     std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
1989
1990     /* Get force params for normal and free energy perturbation
1991      * studies, as determined by types!
1992      */
1993     InteractionOfType param(atomIndices, forceParam, "");
1994
1995     std::vector<InteractionOfType>::iterator foundAParameter = bondtype[ftype].interactionTypes.end();
1996     std::vector<InteractionOfType>::iterator foundBParameter = bondtype[ftype].interactionTypes.end();
1997     if (bBonded)
1998     {
1999         if (NRFPA(ftype) == 0)
2000         {
2001             bFoundA = true;
2002         }
2003         else
2004         {
2005             foundAParameter = defaultInteractionsOfType(ftype, bondtype, atomTypes, &nparam_defA);
2006             if (foundAParameter != bondtype[ftype].interactionTypes.end())
2007             {
2008                 /* Copy the A-state and B-state default parameters. */
2009                 GMX_ASSERT(NRFPA(ftype) + NRFPB(ftype) <= MAXFORCEPARAM,
2010                            "Bonded interactions may have at most 12 parameters");
2011                 gmx::ArrayRef<const real> defaultParam = foundAParameter->forceParam();
2012                 for (int j = 0; (j < NRFPA(ftype) + NRFPB(ftype)); j++)
2013                 {
2014                     param.setForceParameter(j, defaultParam[j]);
2015                 }
2016                 bFoundA = true;
2017             }
2018         }
2019
2020         if (NRFPB(ftype) == 0)
2021         {
2022             bFoundB = true;
2023         }
2024         else
2025         {
2026             foundBParameter = defaultInteractionsOfType(ftype, bondtype, atomTypesB, &nparam_defB);
2027             if (foundBParameter != bondtype[ftype].interactionTypes.end())
2028             {
2029                 /* Copy only the B-state default parameters */
2030                 gmx::ArrayRef<const real> defaultParam = foundBParameter->forceParam();
2031                 for (int j = NRFPA(ftype); (j < NRFP(ftype)); j++)
2032                 {
2033                     param.setForceParameter(j, defaultParam[j]);
2034                 }
2035                 bFoundB = true;
2036             }
2037         }
2038     }
2039     else if (ftype == F_LJ14)
2040     {
2041         bFoundA = default_nb_params(ftype, bondtype, at, &param, 0, FALSE, bGenPairs);
2042         bFoundB = default_nb_params(ftype, bondtype, at, &param, 0, TRUE, bGenPairs);
2043     }
2044     else if (ftype == F_LJC14_Q)
2045     {
2046         /* Fill in the A-state charges as default parameters */
2047         param.setForceParameter(0, fudgeQQ);
2048         param.setForceParameter(1, at->atom[param.ai()].q);
2049         param.setForceParameter(2, at->atom[param.aj()].q);
2050         /* The default LJ parameters are the standard 1-4 parameters */
2051         bFoundA = default_nb_params(F_LJ14, bondtype, at, &param, 3, FALSE, bGenPairs);
2052         bFoundB = TRUE;
2053     }
2054     else if (ftype == F_LJC_PAIRS_NB)
2055     {
2056         /* Defaults are not supported here */
2057         bFoundA = FALSE;
2058         bFoundB = TRUE;
2059     }
2060     else
2061     {
2062         gmx_incons("Unknown function type in push_bond");
2063     }
2064
2065     if (nread > nral_fmt)
2066     {
2067         /* Manually specified parameters - in this case we discard multiple torsion info! */
2068
2069         strcpy(format, asformat[nral_fmt - 1]);
2070         strcat(format, ccformat);
2071
2072         nread = sscanf(line,
2073                        format,
2074                        &cc[0],
2075                        &cc[1],
2076                        &cc[2],
2077                        &cc[3],
2078                        &cc[4],
2079                        &cc[5],
2080                        &cc[6],
2081                        &cc[7],
2082                        &cc[8],
2083                        &cc[9],
2084                        &cc[10],
2085                        &cc[11],
2086                        &cc[12]);
2087
2088         if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
2089         {
2090             /* We only have to issue a warning if these atoms are perturbed! */
2091             bool                     bPert      = false;
2092             gmx::ArrayRef<const int> paramAtoms = param.atoms();
2093             for (int j = 0; (j < nral); j++)
2094             {
2095                 bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]);
2096             }
2097
2098             if (bPert && *bWarn_copy_A_B)
2099             {
2100                 auto message = gmx::formatString(
2101                         "Some parameters for bonded interaction involving "
2102                         "perturbed atoms are specified explicitly in "
2103                         "state A, but not B - copying A to B");
2104                 warning(wi, message);
2105                 *bWarn_copy_A_B = FALSE;
2106             }
2107
2108             /* If only the A parameters were specified, copy them to the B state */
2109             /* The B-state parameters correspond to the first nrfpB
2110              * A-state parameters.
2111              */
2112             for (int j = 0; (j < NRFPB(ftype)); j++)
2113             {
2114                 cc[nread++] = cc[j];
2115             }
2116         }
2117
2118         /* If nread was 0 or EOF, no parameters were read => use defaults.
2119          * If nread was nrfpA we copied above so nread=nrfp.
2120          * If nread was nrfp we are cool.
2121          * For F_LJC14_Q we allow supplying fudgeQQ only.
2122          * Anything else is an error!
2123          */
2124         if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) && !(ftype == F_LJC14_Q && nread == 1))
2125         {
2126             auto message = gmx::formatString(
2127                     "Incorrect number of parameters - found %d, expected %d "
2128                     "or %d for %s (after the function type).",
2129                     nread,
2130                     NRFPA(ftype),
2131                     NRFP(ftype),
2132                     interaction_function[ftype].longname);
2133             warning_error_and_exit(wi, message, FARGS);
2134         }
2135
2136         for (int j = 0; (j < nread); j++)
2137         {
2138             param.setForceParameter(j, cc[j]);
2139         }
2140         /* Check whether we have to use the defaults */
2141         if (nread == NRFP(ftype))
2142         {
2143             bDef = FALSE;
2144         }
2145     }
2146     else
2147     {
2148         nread = 0;
2149     }
2150     /* nread now holds the number of force parameters read! */
2151
2152     if (bDef)
2153     {
2154         /* Use defaults */
2155         /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
2156         if (ftype == F_PDIHS)
2157         {
2158             if ((nparam_defA != nparam_defB)
2159                 || ((nparam_defA > 1 || nparam_defB > 1) && (foundAParameter != foundBParameter)))
2160             {
2161                 auto message = gmx::formatString(
2162                         "Cannot automatically perturb a torsion with multiple terms to different "
2163                         "form.\n"
2164                         "Please specify perturbed parameters manually for this torsion in your "
2165                         "topology!");
2166                 warning_error(wi, message);
2167             }
2168         }
2169
2170         if (nread > 0 && nread < NRFPA(ftype))
2171         {
2172             /* Issue an error, do not use defaults */
2173             auto message = gmx::formatString(
2174                     "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
2175             warning_error(wi, message);
2176         }
2177
2178         if (nread == 0 || nread == EOF)
2179         {
2180             if (!bFoundA)
2181             {
2182                 if (interaction_function[ftype].flags & IF_VSITE)
2183                 {
2184                     for (int j = 0; j < MAXFORCEPARAM; j++)
2185                     {
2186                         param.setForceParameter(j, NOTSET);
2187                     }
2188                     if (bSwapParity)
2189                     {
2190                         /* flag to swap parity of vsi  te construction */
2191                         param.setForceParameter(1, -1);
2192                     }
2193                 }
2194                 else
2195                 {
2196                     if (bZero)
2197                     {
2198                         fprintf(stderr,
2199                                 "NOTE: No default %s types, using zeroes\n",
2200                                 interaction_function[ftype].longname);
2201                     }
2202                     else
2203                     {
2204                         auto message = gmx::formatString("No default %s types",
2205                                                          interaction_function[ftype].longname);
2206                         warning_error(wi, message);
2207                     }
2208                 }
2209             }
2210             else
2211             {
2212                 if (bSwapParity)
2213                 {
2214                     switch (ftype)
2215                     {
2216                         case F_VSITE3FAD: param.setForceParameter(0, 360 - param.c0()); break;
2217                         case F_VSITE3OUT: param.setForceParameter(2, -param.c2()); break;
2218                     }
2219                 }
2220             }
2221             if (!bFoundB)
2222             {
2223                 /* We only have to issue a warning if these atoms are perturbed! */
2224                 bool                     bPert      = false;
2225                 gmx::ArrayRef<const int> paramAtoms = param.atoms();
2226                 for (int j = 0; (j < nral); j++)
2227                 {
2228                     bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]);
2229                 }
2230
2231                 if (bPert)
2232                 {
2233                     auto message = gmx::formatString(
2234                             "No default %s types for perturbed atoms, "
2235                             "using normal values",
2236                             interaction_function[ftype].longname);
2237                     warning(wi, message);
2238                 }
2239             }
2240         }
2241     }
2242
2243     gmx::ArrayRef<const real> paramValue = param.forceParam();
2244     if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ) && paramValue[5] != paramValue[2])
2245     {
2246         auto message = gmx::formatString("%s multiplicity can not be perturbed %f!=%f",
2247                                          interaction_function[ftype].longname,
2248                                          paramValue[2],
2249                                          paramValue[5]);
2250         warning_error_and_exit(wi, message, FARGS);
2251     }
2252
2253     if (IS_TABULATED(ftype) && param.c0() != param.c2())
2254     {
2255         auto message = gmx::formatString("%s table number can not be perturbed %d!=%d",
2256                                          interaction_function[ftype].longname,
2257                                          gmx::roundToInt(param.c0()),
2258                                          gmx::roundToInt(param.c0()));
2259         warning_error_and_exit(wi, message, FARGS);
2260     }
2261
2262     /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2263     if (ftype == F_RBDIHS)
2264     {
2265
2266         int nr = 0;
2267         for (int i = 0; i < NRFP(ftype); i++)
2268         {
2269             if (paramValue[i] != 0.0)
2270             {
2271                 nr++;
2272             }
2273         }
2274         if (nr == 0)
2275         {
2276             return;
2277         }
2278     }
2279
2280     /* Put the values in the appropriate arrays */
2281     add_param_to_list(&bond[ftype], param);
2282
2283     /* Push additional torsions from FF for ftype==9 if we have them.
2284      * We have already checked that the A/B states do not differ in this case,
2285      * so we do not have to double-check that again, or the vsite stuff.
2286      * In addition, those torsions cannot be automatically perturbed.
2287      */
2288     if (bDef && ftype == F_PDIHS)
2289     {
2290         for (int i = 1; i < nparam_defA; i++)
2291         {
2292             /* Advance pointer! */
2293             foundAParameter += 2;
2294             gmx::ArrayRef<const real> forceParam = foundAParameter->forceParam();
2295             for (int j = 0; j < (NRFPA(ftype) + NRFPB(ftype)); j++)
2296             {
2297                 param.setForceParameter(j, forceParam[j]);
2298             }
2299             /* And push the next term for this torsion */
2300             add_param_to_list(&bond[ftype], param);
2301         }
2302     }
2303 }
2304
2305 void push_cmap(Directive                         d,
2306                gmx::ArrayRef<InteractionsOfType> bondtype,
2307                gmx::ArrayRef<InteractionsOfType> bond,
2308                t_atoms*                          at,
2309                PreprocessingAtomTypes*           atypes,
2310                char*                             line,
2311                warninp*                          wi)
2312 {
2313     const char* aaformat[MAXATOMLIST + 1] = {
2314         "%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"
2315     };
2316
2317     int  ftype, nral, nread, ncmap_params;
2318     int  cmap_type;
2319     int  aa[MAXATOMLIST + 1];
2320     bool bFound;
2321
2322     ftype        = ifunc_index(d, 1);
2323     nral         = NRAL(ftype);
2324     ncmap_params = 0;
2325
2326     nread = sscanf(line, aaformat[nral - 1], &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2327
2328     if (nread < nral)
2329     {
2330         too_few(wi);
2331         return;
2332     }
2333     else if (nread == nral)
2334     {
2335         ftype = ifunc_index(d, 1);
2336     }
2337
2338     /* Check for double atoms and atoms out of bounds */
2339     for (int i = 0; i < nral; i++)
2340     {
2341         if (aa[i] < 1 || aa[i] > at->nr)
2342         {
2343             auto message = gmx::formatString(
2344                     "Atom index (%d) in %s out of bounds (1-%d).\n"
2345                     "This probably means that you have inserted topology section \"%s\"\n"
2346                     "in a part belonging to a different molecule than you intended to.\n"
2347                     "In that case move the \"%s\" section to the right molecule.",
2348                     aa[i],
2349                     enumValueToString(d),
2350                     at->nr,
2351                     enumValueToString(d),
2352                     enumValueToString(d));
2353             warning_error_and_exit(wi, message, FARGS);
2354         }
2355
2356         for (int j = i + 1; (j < nral); j++)
2357         {
2358             if (aa[i] == aa[j])
2359             {
2360                 auto message = gmx::formatString(
2361                         "Duplicate atom index (%d) in %s", aa[i], enumValueToString(d));
2362                 warning_error(wi, message);
2363             }
2364         }
2365     }
2366
2367     /* default force parameters  */
2368     std::vector<int> atoms;
2369     for (int j = 0; (j < nral); j++)
2370     {
2371         atoms.emplace_back(aa[j] - 1);
2372     }
2373     std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2374     InteractionOfType               param(atoms, forceParam, "");
2375     /* Get the cmap type for this cmap angle */
2376     bFound = default_cmap_params(bondtype, at, atypes, &param, FALSE, &cmap_type, &ncmap_params, wi);
2377
2378     /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2379     if (bFound && ncmap_params == 1)
2380     {
2381         /* Put the values in the appropriate arrays */
2382         param.setForceParameter(0, cmap_type);
2383         add_param_to_list(&bond[ftype], param);
2384     }
2385     else
2386     {
2387         /* This is essentially the same check as in default_cmap_params() done one more time */
2388         auto message =
2389                 gmx::formatString("Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2390                                   param.ai() + 1,
2391                                   param.aj() + 1,
2392                                   param.ak() + 1,
2393                                   param.al() + 1,
2394                                   param.am() + 1);
2395         warning_error_and_exit(wi, message, FARGS);
2396     }
2397 }
2398
2399
2400 void push_vsitesn(Directive d, gmx::ArrayRef<InteractionsOfType> bond, t_atoms* at, char* line, warninp* wi)
2401 {
2402     char*   ptr;
2403     int     type, ftype, n, ret, nj, a;
2404     int*    atc    = nullptr;
2405     double *weight = nullptr, weight_tot;
2406
2407     std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2408     ptr                                        = line;
2409     ret                                        = sscanf(ptr, "%d%n", &a, &n);
2410     ptr += n;
2411     if (ret == 0)
2412     {
2413         auto message =
2414                 gmx::formatString("Expected an atom index in section \"%s\"", enumValueToString(d));
2415         warning_error_and_exit(wi, message, FARGS);
2416     }
2417
2418     sscanf(ptr, "%d%n", &type, &n);
2419     ptr += n;
2420     ftype         = ifunc_index(d, type);
2421     int firstAtom = a - 1;
2422
2423     weight_tot = 0;
2424     nj         = 0;
2425     do
2426     {
2427         ret = sscanf(ptr, "%d%n", &a, &n);
2428         ptr += n;
2429         if (ret > 0)
2430         {
2431             if (nj % 20 == 0)
2432             {
2433                 srenew(atc, nj + 20);
2434                 srenew(weight, nj + 20);
2435             }
2436             atc[nj] = a - 1;
2437             switch (type)
2438             {
2439                 case 1: weight[nj] = 1; break;
2440                 case 2:
2441                     /* Here we use the A-state mass as a parameter.
2442                      * Note that the B-state mass has no influence.
2443                      */
2444                     weight[nj] = at->atom[atc[nj]].m;
2445                     break;
2446                 case 3:
2447                     weight[nj] = -1;
2448                     ret        = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2449                     ptr += n;
2450                     if (weight[nj] < 0)
2451                     {
2452                         auto message = gmx::formatString(
2453                                 "No weight or negative weight found for vsiten "
2454                                 "constructing atom %d (atom index %d)",
2455                                 nj + 1,
2456                                 atc[nj] + 1);
2457                         warning_error_and_exit(wi, message, FARGS);
2458                     }
2459                     break;
2460                 default:
2461                     auto message = gmx::formatString("Unknown vsiten type %d", type);
2462                     warning_error_and_exit(wi, message, FARGS);
2463             }
2464             weight_tot += weight[nj];
2465             nj++;
2466         }
2467     } while (ret > 0);
2468
2469     if (nj == 0)
2470     {
2471         auto message = gmx::formatString("Expected more than one atom index in section \"%s\"",
2472                                          enumValueToString(d));
2473         warning_error_and_exit(wi, message, FARGS);
2474     }
2475
2476     if (weight_tot == 0)
2477     {
2478         warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
2479     }
2480
2481     for (int j = 0; j < nj; j++)
2482     {
2483         std::vector<int> atoms = { firstAtom, atc[j] };
2484         forceParam[0]          = nj;
2485         forceParam[1]          = weight[j] / weight_tot;
2486         /* Put the values in the appropriate arrays */
2487         add_param_to_list(&bond[ftype], InteractionOfType(atoms, forceParam));
2488     }
2489
2490     sfree(atc);
2491     sfree(weight);
2492 }
2493
2494 void push_mol(gmx::ArrayRef<MoleculeInformation> mols, char* pline, int* whichmol, int* nrcopies, warninp* wi)
2495 {
2496     char type[STRLEN];
2497
2498     if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2499     {
2500         too_few(wi);
2501         return;
2502     }
2503
2504     /* Search moleculename.
2505      * Here we originally only did case insensitive matching. But because
2506      * some PDB files can have many chains and use case to generate more
2507      * chain-identifiers, which in turn end up in our moleculetype name,
2508      * we added support for case-sensitivity.
2509      */
2510     int nrcs    = 0;
2511     int nrci    = 0;
2512     int matchci = -1;
2513     int matchcs = -1;
2514     int i       = 0;
2515     for (const auto& mol : mols)
2516     {
2517         if (strcmp(type, *(mol.name)) == 0)
2518         {
2519             nrcs++;
2520             matchcs = i;
2521         }
2522         if (gmx_strcasecmp(type, *(mol.name)) == 0)
2523         {
2524             nrci++;
2525             matchci = i;
2526         }
2527         i++;
2528     }
2529
2530     if (nrcs == 1)
2531     {
2532         // select the case sensitive match
2533         *whichmol = matchcs;
2534     }
2535     else
2536     {
2537         // avoid matching case-insensitive when we have multiple matches
2538         if (nrci > 1)
2539         {
2540             auto message = gmx::formatString(
2541                     "For moleculetype '%s' in [ system ] %d case insensitive "
2542                     "matches, but %d case sensitive matches were found. Check "
2543                     "the case of the characters in the moleculetypes.",
2544                     type,
2545                     nrci,
2546                     nrcs);
2547             warning_error_and_exit(wi, message, FARGS);
2548         }
2549         if (nrci == 1)
2550         {
2551             // select the unique case insensitive match
2552             *whichmol = matchci;
2553         }
2554         else
2555         {
2556             auto message = gmx::formatString("No such moleculetype %s", type);
2557             warning_error_and_exit(wi, message, FARGS);
2558         }
2559     }
2560 }
2561
2562 void push_excl(char* line, gmx::ArrayRef<gmx::ExclusionBlock> b2, warninp* wi)
2563 {
2564     int  i, j;
2565     int  n;
2566     char base[STRLEN], format[STRLEN];
2567
2568     if (sscanf(line, "%d", &i) == 0)
2569     {
2570         return;
2571     }
2572
2573     if ((1 <= i) && (i <= b2.ssize()))
2574     {
2575         i--;
2576     }
2577     else
2578     {
2579         return;
2580     }
2581     strcpy(base, "%*d");
2582     do
2583     {
2584         strcpy(format, base);
2585         strcat(format, "%d");
2586         n = sscanf(line, format, &j);
2587         if (n == 1)
2588         {
2589             if ((1 <= j) && (j <= b2.ssize()))
2590             {
2591                 j--;
2592                 b2[i].atomNumber.push_back(j);
2593                 /* also add the reverse exclusion! */
2594                 b2[j].atomNumber.push_back(i);
2595                 strcat(base, "%*d");
2596             }
2597             else
2598             {
2599                 auto message = gmx::formatString("Invalid Atomnr j: %d, b2->nr: %zu\n", j, b2.size());
2600                 warning_error_and_exit(wi, message, FARGS);
2601             }
2602         }
2603     } while (n == 1);
2604 }
2605
2606 int add_atomtype_decoupled(t_symtab* symtab, PreprocessingAtomTypes* at, t_nbparam*** nbparam, t_nbparam*** pair)
2607 {
2608     t_atom atom;
2609     int    nr;
2610
2611     /* Define an atom type with all parameters set to zero (no interactions) */
2612     atom.q = 0.0;
2613     atom.m = 0.0;
2614     /* Type for decoupled atoms could be anything,
2615      * this should be changed automatically later when required.
2616      */
2617     atom.ptype = ParticleType::Atom;
2618
2619     std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2620     nr = at->addType(symtab, atom, "decoupled", InteractionOfType({}, forceParam, ""), -1, 0);
2621
2622     /* Add space in the non-bonded parameters matrix */
2623     realloc_nb_params(at, nbparam, pair);
2624
2625     return nr;
2626 }
2627
2628 static void convert_pairs_to_pairsQ(gmx::ArrayRef<InteractionsOfType> interactions, real fudgeQQ, t_atoms* atoms)
2629 {
2630     /* Add the pair list to the pairQ list */
2631     std::vector<InteractionOfType> paramnew;
2632
2633     gmx::ArrayRef<const InteractionOfType> paramp1 = interactions[F_LJ14].interactionTypes;
2634     gmx::ArrayRef<const InteractionOfType> paramp2 = interactions[F_LJC14_Q].interactionTypes;
2635
2636     /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2637        it may be possible to just ADD the converted F_LJ14 array
2638        to the old F_LJC14_Q array, but since we have to create
2639        a new sized memory structure, better just to deep copy it all.
2640      */
2641
2642
2643     for (const auto& param : paramp2)
2644     {
2645         paramnew.emplace_back(param);
2646     }
2647
2648     for (const auto& param : paramp1)
2649     {
2650         std::vector<real> forceParam = {
2651             fudgeQQ, atoms->atom[param.ai()].q, atoms->atom[param.aj()].q, param.c0(), param.c1()
2652         };
2653         paramnew.emplace_back(param.atoms(), forceParam, "");
2654     }
2655
2656     /* now assign the new data to the F_LJC14_Q structure */
2657     interactions[F_LJC14_Q].interactionTypes = paramnew;
2658
2659     /* Empty the LJ14 pairlist */
2660     interactions[F_LJ14].interactionTypes.clear();
2661 }
2662
2663 static void generate_LJCpairsNB(MoleculeInformation* mol, int nb_funct, InteractionsOfType* nbp, warninp* wi)
2664 {
2665     int     n, ntype;
2666     t_atom* atom;
2667
2668     n    = mol->atoms.nr;
2669     atom = mol->atoms.atom;
2670
2671     ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->size())));
2672     GMX_ASSERT(ntype * ntype == gmx::ssize(*nbp),
2673                "Number of pairs of generated non-bonded parameters should be a perfect square");
2674
2675     /* Add a pair interaction for all non-excluded atom pairs */
2676     const auto& excls = mol->excls;
2677     for (int i = 0; i < n; i++)
2678     {
2679         for (int j = i + 1; j < n; j++)
2680         {
2681             bool pairIsExcluded = false;
2682             for (const int atomK : excls[i])
2683             {
2684                 if (atomK == j)
2685                 {
2686                     pairIsExcluded = true;
2687                 }
2688             }
2689             if (!pairIsExcluded)
2690             {
2691                 if (nb_funct != F_LJ)
2692                 {
2693                     auto message = gmx::formatString(
2694                             "Can only generate non-bonded pair interactions "
2695                             "for Van der Waals type Lennard-Jones");
2696                     warning_error_and_exit(wi, message, FARGS);
2697                 }
2698                 std::vector<int>  atoms      = { i, j };
2699                 std::vector<real> forceParam = {
2700                     atom[i].q,
2701                     atom[j].q,
2702                     nbp->interactionTypes[ntype * atom[i].type + atom[j].type].c0(),
2703                     nbp->interactionTypes[ntype * atom[i].type + atom[j].type].c1()
2704                 };
2705                 add_param_to_list(&mol->interactions[F_LJC_PAIRS_NB], InteractionOfType(atoms, forceParam));
2706             }
2707         }
2708     }
2709 }
2710
2711 static void set_excl_all(gmx::ListOfLists<int>* excl)
2712 {
2713     /* Get rid of the current exclusions and exclude all atom pairs */
2714     const int        numAtoms = excl->ssize();
2715     std::vector<int> exclusionsForAtom(numAtoms);
2716     for (int i = 0; i < numAtoms; i++)
2717     {
2718         exclusionsForAtom[i] = i;
2719     }
2720     excl->clear();
2721     for (int i = 0; i < numAtoms; i++)
2722     {
2723         excl->pushBack(exclusionsForAtom);
2724     }
2725 }
2726
2727 static void decouple_atoms(t_atoms*    atoms,
2728                            int         atomtype_decouple,
2729                            int         couple_lam0,
2730                            int         couple_lam1,
2731                            const char* mol_name,
2732                            warninp*    wi)
2733 {
2734     int i;
2735
2736     for (i = 0; i < atoms->nr; i++)
2737     {
2738         t_atom* atom;
2739
2740         atom = &atoms->atom[i];
2741
2742         if (atom->qB != atom->q || atom->typeB != atom->type)
2743         {
2744             auto message = gmx::formatString(
2745                     "Atom %d in molecule type '%s' has different A and B state "
2746                     "charges and/or atom types set in the topology file as well "
2747                     "as through the mdp option '%s'. You can not use both "
2748                     "these methods simultaneously.",
2749                     i + 1,
2750                     mol_name,
2751                     "couple-moltype");
2752             warning_error_and_exit(wi, message, FARGS);
2753         }
2754
2755         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2756         {
2757             atom->q = 0.0;
2758         }
2759         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2760         {
2761             atom->type = atomtype_decouple;
2762         }
2763         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2764         {
2765             atom->qB = 0.0;
2766         }
2767         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2768         {
2769             atom->typeB = atomtype_decouple;
2770         }
2771     }
2772 }
2773
2774 void convert_moltype_couple(MoleculeInformation* mol,
2775                             int                  atomtype_decouple,
2776                             real                 fudgeQQ,
2777                             int                  couple_lam0,
2778                             int                  couple_lam1,
2779                             bool                 bCoupleIntra,
2780                             int                  nb_funct,
2781                             InteractionsOfType*  nbp,
2782                             warninp*             wi)
2783 {
2784     convert_pairs_to_pairsQ(mol->interactions, fudgeQQ, &mol->atoms);
2785     if (!bCoupleIntra)
2786     {
2787         generate_LJCpairsNB(mol, nb_funct, nbp, wi);
2788         set_excl_all(&mol->excls);
2789     }
2790     decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1, *mol->name, wi);
2791 }