src/gromacs/gmxpreprocess/toppush.cpp

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