src/gromacs/gmxpreprocess/toppush.cpp

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