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