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