src/gromacs/gmxpreprocess/toppush.cpp

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