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