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