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