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