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