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