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