src/gromacs/gmxpreprocess/toppush.cpp

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