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