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