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