src/gromacs/gmxpreprocess/toppush.cpp

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