src/gromacs/gmxpreprocess/toppush.cpp

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