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 static 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             GMX_ASSERT(j < MAXATOMLIST + 1, "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1886             if (aa[i] == aa[j])
1887             {
1888                 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1889                 if (ftype == F_ANGRES)
1890                 {
1891                     /* Since the angle restraints uses 2 pairs of atoms to
1892                      * defines an angle between vectors, it can be useful
1893                      * to use one atom twice, so we only issue a note here.
1894                      */
1895                     warning_note(wi, errbuf);
1896                 }
1897                 else
1898                 {
1899                     warning_error(wi, errbuf);
1900                 }
1901             }
1902         }
1903     }
1904
1905     /* default force parameters  */
1906     for (j = 0; (j < MAXATOMLIST); j++)
1907     {
1908         param.a[j] = aa[j]-1;
1909     }
1910     for (j = 0; (j < MAXFORCEPARAM); j++)
1911     {
1912         param.c[j] = 0.0;
1913     }
1914
1915     /* Get force params for normal and free energy perturbation
1916      * studies, as determined by types!
1917      */
1918
1919     if (bBonded)
1920     {
1921         bFoundA = default_params(ftype, bondtype, at, atype, &param, FALSE, &param_defA, &nparam_defA);
1922         if (bFoundA)
1923         {
1924             /* Copy the A-state and B-state default parameters. */
1925             GMX_ASSERT(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM, "Bonded interactions may have at most 12 parameters");
1926             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
1927             {
1928                 param.c[j] = param_defA->c[j];
1929             }
1930         }
1931         bFoundB = default_params(ftype, bondtype, at, atype, &param, TRUE, &param_defB, &nparam_defB);
1932         if (bFoundB)
1933         {
1934             /* Copy only the B-state default parameters */
1935             for (j = NRFPA(ftype); (j < NRFP(ftype)); j++)
1936             {
1937                 param.c[j] = param_defB->c[j];
1938             }
1939         }
1940     }
1941     else if (ftype == F_LJ14)
1942     {
1943         bFoundA = default_nb_params(ftype, bondtype, at, &param, 0, FALSE, bGenPairs);
1944         bFoundB = default_nb_params(ftype, bondtype, at, &param, 0, TRUE, bGenPairs);
1945     }
1946     else if (ftype == F_LJC14_Q)
1947     {
1948         param.c[0] = fudgeQQ;
1949         /* Fill in the A-state charges as default parameters */
1950         param.c[1] = at->atom[param.a[0]].q;
1951         param.c[2] = at->atom[param.a[1]].q;
1952         /* The default LJ parameters are the standard 1-4 parameters */
1953         bFoundA = default_nb_params(F_LJ14, bondtype, at, &param, 3, FALSE, bGenPairs);
1954         bFoundB = TRUE;
1955     }
1956     else if (ftype == F_LJC_PAIRS_NB)
1957     {
1958         /* Defaults are not supported here */
1959         bFoundA = FALSE;
1960         bFoundB = TRUE;
1961     }
1962     else
1963     {
1964         gmx_incons("Unknown function type in push_bond");
1965     }
1966
1967     if (nread > nral_fmt)
1968     {
1969         /* Manually specified parameters - in this case we discard multiple torsion info! */
1970
1971         strcpy(format, asformat[nral_fmt-1]);
1972         strcat(format, ccformat);
1973
1974         nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5],
1975                        &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
1976
1977         if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1978         {
1979             /* We only have to issue a warning if these atoms are perturbed! */
1980             bPert = FALSE;
1981             for (j = 0; (j < nral); j++)
1982             {
1983                 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1984             }
1985
1986             if (bPert && *bWarn_copy_A_B)
1987             {
1988                 sprintf(errbuf,
1989                         "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1990                 warning(wi, errbuf);
1991                 *bWarn_copy_A_B = FALSE;
1992             }
1993
1994             /* If only the A parameters were specified, copy them to the B state */
1995             /* The B-state parameters correspond to the first nrfpB
1996              * A-state parameters.
1997              */
1998             for (j = 0; (j < NRFPB(ftype)); j++)
1999             {
2000                 cc[nread++] = cc[j];
2001             }
2002         }
2003
2004         /* If nread was 0 or EOF, no parameters were read => use defaults.
2005          * If nread was nrfpA we copied above so nread=nrfp.
2006          * If nread was nrfp we are cool.
2007          * For F_LJC14_Q we allow supplying fudgeQQ only.
2008          * Anything else is an error!
2009          */
2010         if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
2011             !(ftype == F_LJC14_Q && nread == 1))
2012         {
2013             sprintf(errbuf, "Incorrect number of parameters - found %d, expected %d or %d for %s.",
2014                     nread, NRFPA(ftype), NRFP(ftype), interaction_function[ftype].longname);
2015             warning_error_and_exit(wi, errbuf, FARGS);
2016         }
2017
2018         for (j = 0; (j < nread); j++)
2019         {
2020             param.c[j] = cc[j];
2021         }
2022
2023         /* Check whether we have to use the defaults */
2024         if (nread == NRFP(ftype))
2025         {
2026             bDef = FALSE;
2027         }
2028     }
2029     else
2030     {
2031         nread = 0;
2032     }
2033     /* nread now holds the number of force parameters read! */
2034
2035     if (bDef)
2036     {
2037         /* Use defaults */
2038         /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
2039         if (ftype == F_PDIHS)
2040         {
2041             if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB)))
2042             {
2043                 sprintf(errbuf,
2044                         "Cannot automatically perturb a torsion with multiple terms to different form.\n"
2045                         "Please specify perturbed parameters manually for this torsion in your topology!");
2046                 warning_error(wi, errbuf);
2047             }
2048         }
2049
2050         if (nread > 0 && nread < NRFPA(ftype))
2051         {
2052             /* Issue an error, do not use defaults */
2053             sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
2054             warning_error(wi, errbuf);
2055         }
2056
2057         if (nread == 0 || nread == EOF)
2058         {
2059             if (!bFoundA)
2060             {
2061                 if (interaction_function[ftype].flags & IF_VSITE)
2062                 {
2063                     /* set them to NOTSET, will be calculated later */
2064                     for (j = 0; (j < MAXFORCEPARAM); j++)
2065                     {
2066                         param.c[j] = NOTSET;
2067                     }
2068
2069                     if (bSwapParity)
2070                     {
2071                         param.c1() = -1; /* flag to swap parity of vsite construction */
2072                     }
2073                 }
2074                 else
2075                 {
2076                     if (bZero)
2077                     {
2078                         fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
2079                                 interaction_function[ftype].longname);
2080                     }
2081                     else
2082                     {
2083                         sprintf(errbuf, "No default %s types", interaction_function[ftype].longname);
2084                         warning_error(wi, errbuf);
2085                     }
2086                 }
2087             }
2088             else
2089             {
2090                 if (bSwapParity)
2091                 {
2092                     switch (ftype)
2093                     {
2094                         case F_VSITE3FAD:
2095                             param.c0() = 360 - param.c0();
2096                             break;
2097                         case F_VSITE3OUT:
2098                             param.c2() = -param.c2();
2099                             break;
2100                     }
2101                 }
2102             }
2103             if (!bFoundB)
2104             {
2105                 /* We only have to issue a warning if these atoms are perturbed! */
2106                 bPert = FALSE;
2107                 for (j = 0; (j < nral); j++)
2108                 {
2109                     bPert = bPert || PERTURBED(at->atom[param.a[j]]);
2110                 }
2111
2112                 if (bPert)
2113                 {
2114                     sprintf(errbuf, "No default %s types for perturbed atoms, "
2115                             "using normal values", interaction_function[ftype].longname);
2116                     warning(wi, errbuf);
2117                 }
2118             }
2119         }
2120     }
2121
2122     if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
2123         && param.c[5] != param.c[2])
2124     {
2125         sprintf(errbuf, "%s multiplicity can not be perturbed %f!=%f",
2126                 interaction_function[ftype].longname,
2127                 param.c[2], param.c[5]);
2128         warning_error_and_exit(wi, errbuf, FARGS);
2129     }
2130
2131     if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
2132     {
2133         sprintf(errbuf, "%s table number can not be perturbed %d!=%d",
2134                 interaction_function[ftype].longname,
2135                 (int)(param.c[0]+0.5), (int)(param.c[2]+0.5));
2136         warning_error_and_exit(wi, errbuf, FARGS);
2137     }
2138
2139     /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2140     if (ftype == F_RBDIHS)
2141     {
2142         nr = 0;
2143         for (i = 0; i < NRFP(ftype); i++)
2144         {
2145             if (param.c[i] != 0)
2146             {
2147                 nr++;
2148             }
2149         }
2150         if (nr == 0)
2151         {
2152             return;
2153         }
2154     }
2155
2156     /* Put the values in the appropriate arrays */
2157     add_param_to_list (&bond[ftype], &param);
2158
2159     /* Push additional torsions from FF for ftype==9 if we have them.
2160      * We have already checked that the A/B states do not differ in this case,
2161      * so we do not have to double-check that again, or the vsite stuff.
2162      * In addition, those torsions cannot be automatically perturbed.
2163      */
2164     if (bDef && ftype == F_PDIHS)
2165     {
2166         for (i = 1; i < nparam_defA; i++)
2167         {
2168             /* Advance pointer! */
2169             param_defA += 2;
2170             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
2171             {
2172                 param.c[j] = param_defA->c[j];
2173             }
2174             /* And push the next term for this torsion */
2175             add_param_to_list (&bond[ftype], &param);
2176         }
2177     }
2178 }
2179
2180 void push_cmap(directive d, t_params bondtype[], t_params bond[],
2181                t_atoms *at, gpp_atomtype_t atype, char *line,
2182                warninp_t wi)
2183 {
2184     const char *aaformat[MAXATOMLIST+1] =
2185     {
2186         "%d",
2187         "%d%d",
2188         "%d%d%d",
2189         "%d%d%d%d",
2190         "%d%d%d%d%d",
2191         "%d%d%d%d%d%d",
2192         "%d%d%d%d%d%d%d"
2193     };
2194
2195     int      i, j, ftype, nral, nread, ncmap_params;
2196     int      cmap_type;
2197     int      aa[MAXATOMLIST+1];
2198     char     errbuf[STRLEN];
2199     gmx_bool bFound;
2200     t_param  param;
2201
2202     ftype        = ifunc_index(d, 1);
2203     nral         = NRAL(ftype);
2204     ncmap_params = 0;
2205
2206     nread = sscanf(line, aaformat[nral-1],
2207                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2208
2209     if (nread < nral)
2210     {
2211         too_few(wi);
2212         return;
2213     }
2214     else if (nread == nral)
2215     {
2216         ftype = ifunc_index(d, 1);
2217     }
2218
2219     /* Check for double atoms and atoms out of bounds */
2220     for (i = 0; i < nral; i++)
2221     {
2222         if (aa[i] < 1 || aa[i] > at->nr)
2223         {
2224             sprintf(errbuf, "Atom index (%d) in %s out of bounds (1-%d).\n"
2225                     "This probably means that you have inserted topology section \"%s\"\n"
2226                     "in a part belonging to a different molecule than you intended to.\n"
2227                     "In that case move the \"%s\" section to the right molecule.",
2228                     aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2229             warning_error_and_exit(wi, errbuf, FARGS);
2230         }
2231
2232         for (j = i+1; (j < nral); j++)
2233         {
2234             if (aa[i] == aa[j])
2235             {
2236                 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2237                 warning_error(wi, errbuf);
2238             }
2239         }
2240     }
2241
2242     /* default force parameters  */
2243     for (j = 0; (j < MAXATOMLIST); j++)
2244     {
2245         param.a[j] = aa[j]-1;
2246     }
2247     for (j = 0; (j < MAXFORCEPARAM); j++)
2248     {
2249         param.c[j] = 0.0;
2250     }
2251
2252     /* Get the cmap type for this cmap angle */
2253     bFound = default_cmap_params(bondtype, at, atype, &param, FALSE, &cmap_type, &ncmap_params, wi);
2254
2255     /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2256     if (bFound && ncmap_params == 1)
2257     {
2258         /* Put the values in the appropriate arrays */
2259         param.c[0] = cmap_type;
2260         add_param_to_list(&bond[ftype], &param);
2261     }
2262     else
2263     {
2264         /* This is essentially the same check as in default_cmap_params() done one more time */
2265         sprintf(errbuf, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2266                 param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2267         warning_error_and_exit(wi, errbuf, FARGS);
2268     }
2269 }
2270
2271
2272
2273 void push_vsitesn(directive d, t_params bond[],
2274                   t_atoms *at, char *line,
2275                   warninp_t wi)
2276 {
2277     char   *ptr;
2278     int     type, ftype, j, n, ret, nj, a;
2279     int    *atc    = nullptr;
2280     double *weight = nullptr, weight_tot;
2281     t_param param;
2282     char    errbuf[STRLEN];
2283
2284     /* default force parameters  */
2285     for (j = 0; (j < MAXATOMLIST); j++)
2286     {
2287         param.a[j] = NOTSET;
2288     }
2289     for (j = 0; (j < MAXFORCEPARAM); j++)
2290     {
2291         param.c[j] = 0.0;
2292     }
2293
2294     ptr  = line;
2295     ret  = sscanf(ptr, "%d%n", &a, &n);
2296     ptr += n;
2297     if (ret == 0)
2298     {
2299         sprintf(errbuf, "Expected an atom index in section \"%s\"",
2300                 dir2str(d));
2301         warning_error_and_exit(wi, errbuf, FARGS);
2302     }
2303
2304     param.a[0] = a - 1;
2305
2306     sscanf(ptr, "%d%n", &type, &n);
2307     ptr  += n;
2308     ftype = ifunc_index(d, type);
2309
2310     weight_tot = 0;
2311     nj         = 0;
2312     do
2313     {
2314         ret  = sscanf(ptr, "%d%n", &a, &n);
2315         ptr += n;
2316         if (ret > 0)
2317         {
2318             if (nj % 20 == 0)
2319             {
2320                 srenew(atc, nj+20);
2321                 srenew(weight, nj+20);
2322             }
2323             atc[nj] = a - 1;
2324             switch (type)
2325             {
2326                 case 1:
2327                     weight[nj] = 1;
2328                     break;
2329                 case 2:
2330                     /* Here we use the A-state mass as a parameter.
2331                      * Note that the B-state mass has no influence.
2332                      */
2333                     weight[nj] = at->atom[atc[nj]].m;
2334                     break;
2335                 case 3:
2336                     weight[nj] = -1;
2337                     ret        = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2338                     ptr       += n;
2339                     if (weight[nj] < 0)
2340                     {
2341                         sprintf(errbuf, "No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
2342                                 nj+1, atc[nj]+1);
2343                         warning_error_and_exit(wi, errbuf, FARGS);
2344                     }
2345                     break;
2346                 default:
2347                     sprintf(errbuf, "Unknown vsiten type %d", type);
2348                     warning_error_and_exit(wi, errbuf, FARGS);
2349             }
2350             weight_tot += weight[nj];
2351             nj++;
2352         }
2353     }
2354     while (ret > 0);
2355
2356     if (nj == 0)
2357     {
2358         sprintf(errbuf, "Expected more than one atom index in section \"%s\"",
2359                 dir2str(d));
2360         warning_error_and_exit(wi, errbuf, FARGS);
2361     }
2362
2363     if (weight_tot == 0)
2364     {
2365         warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
2366     }
2367
2368     for (j = 0; j < nj; j++)
2369     {
2370         param.a[1] = atc[j];
2371         param.c[0] = nj;
2372         param.c[1] = weight[j]/weight_tot;
2373         /* Put the values in the appropriate arrays */
2374         add_param_to_list (&bond[ftype], &param);
2375     }
2376
2377     sfree(atc);
2378     sfree(weight);
2379 }
2380
2381 void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
2382               int *nrcopies,
2383               warninp_t wi)
2384 {
2385     char type[STRLEN];
2386     char errbuf[STRLEN];
2387
2388     if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2389     {
2390         too_few(wi);
2391         return;
2392     }
2393
2394     /* Search moleculename.
2395      * Here we originally only did case insensitive matching. But because
2396      * some PDB files can have many chains and use case to generate more
2397      * chain-identifiers, which in turn end up in our moleculetype name,
2398      * we added support for case-sensitivity.
2399      */
2400     int nrcs    = 0;
2401     int nrci    = 0;
2402     int matchci = -1;
2403     int matchcs = -1;
2404     for (int i = 0; i < nrmols; i++)
2405     {
2406         if (strcmp(type, *(mols[i].name)) == 0)
2407         {
2408             nrcs++;
2409             matchcs = i;
2410         }
2411         if (gmx_strcasecmp(type, *(mols[i].name)) == 0)
2412         {
2413             nrci++;
2414             matchci = i;
2415         }
2416     }
2417
2418     if (nrcs == 1)
2419     {
2420         // select the case sensitive match
2421         *whichmol = matchcs;
2422     }
2423     else
2424     {
2425         // avoid matching case-insensitive when we have multiple matches
2426         if (nrci > 1)
2427         {
2428             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);
2429             warning_error_and_exit(wi, errbuf, FARGS);
2430         }
2431         if (nrci == 1)
2432         {
2433             // select the unique case insensitive match
2434             *whichmol = matchci;
2435         }
2436         else
2437         {
2438             sprintf(errbuf, "No such moleculetype %s", type);
2439             warning_error_and_exit(wi, errbuf, FARGS);
2440         }
2441     }
2442 }
2443
2444 void init_block2(t_block2 *b2, int natom)
2445 {
2446     int i;
2447
2448     b2->nr = natom;
2449     snew(b2->nra, b2->nr);
2450     snew(b2->a, b2->nr);
2451     for (i = 0; (i < b2->nr); i++)
2452     {
2453         b2->a[i] = nullptr;
2454     }
2455 }
2456
2457 void done_block2(t_block2 *b2)
2458 {
2459     int i;
2460
2461     if (b2->nr)
2462     {
2463         for (i = 0; (i < b2->nr); i++)
2464         {
2465             sfree(b2->a[i]);
2466         }
2467         sfree(b2->a);
2468         sfree(b2->nra);
2469         b2->nr = 0;
2470     }
2471 }
2472
2473 void push_excl(char *line, t_block2 *b2, warninp_t wi)
2474 {
2475     int  i, j;
2476     int  n;
2477     char base[STRLEN], format[STRLEN];
2478     char errbuf[STRLEN];
2479
2480     if (sscanf(line, "%d", &i) == 0)
2481     {
2482         return;
2483     }
2484
2485     if ((1 <= i) && (i <= b2->nr))
2486     {
2487         i--;
2488     }
2489     else
2490     {
2491         if (debug)
2492         {
2493             fprintf(debug, "Unbound atom %d\n", i-1);
2494         }
2495         return;
2496     }
2497     strcpy(base, "%*d");
2498     do
2499     {
2500         strcpy(format, base);
2501         strcat(format, "%d");
2502         n = sscanf(line, format, &j);
2503         if (n == 1)
2504         {
2505             if ((1 <= j) && (j <= b2->nr))
2506             {
2507                 j--;
2508                 srenew(b2->a[i], ++(b2->nra[i]));
2509                 b2->a[i][b2->nra[i]-1] = j;
2510                 /* also add the reverse exclusion! */
2511                 srenew(b2->a[j], ++(b2->nra[j]));
2512                 b2->a[j][b2->nra[j]-1] = i;
2513                 strcat(base, "%*d");
2514             }
2515             else
2516             {
2517                 sprintf(errbuf, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
2518                 warning_error_and_exit(wi, errbuf, FARGS);
2519             }
2520         }
2521     }
2522     while (n == 1);
2523 }
2524
2525 void b_to_b2(t_blocka *b, t_block2 *b2)
2526 {
2527     int     i;
2528     int     j, a;
2529
2530     for (i = 0; (i < b->nr); i++)
2531     {
2532         for (j = b->index[i]; (j < b->index[i+1]); j++)
2533         {
2534             a = b->a[j];
2535             srenew(b2->a[i], ++b2->nra[i]);
2536             b2->a[i][b2->nra[i]-1] = a;
2537         }
2538     }
2539 }
2540
2541 void b2_to_b(t_block2 *b2, t_blocka *b)
2542 {
2543     int     i, nra;
2544     int     j;
2545
2546     nra = 0;
2547     for (i = 0; (i < b2->nr); i++)
2548     {
2549         b->index[i] = nra;
2550         for (j = 0; (j < b2->nra[i]); j++)
2551         {
2552             b->a[nra+j] = b2->a[i][j];
2553         }
2554         nra += b2->nra[i];
2555     }
2556     /* terminate list */
2557     b->index[i] = nra;
2558 }
2559
2560 static int icomp(const void *v1, const void *v2)
2561 {
2562     return (*((int *) v1))-(*((int *) v2));
2563 }
2564
2565 void merge_excl(t_blocka *excl, t_block2 *b2, warninp_t wi)
2566 {
2567     int     i, k;
2568     int     j;
2569     int     nra;
2570     char    errbuf[STRLEN];
2571
2572     if (!b2->nr)
2573     {
2574         return;
2575     }
2576     else if (b2->nr != excl->nr)
2577     {
2578         sprintf(errbuf, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2579                 b2->nr, excl->nr);
2580         warning_error_and_exit(wi, errbuf, FARGS);
2581     }
2582     else if (debug)
2583     {
2584         fprintf(debug, "Entering merge_excl\n");
2585     }
2586
2587     /* First copy all entries from excl to b2 */
2588     b_to_b2(excl, b2);
2589
2590     /* Count and sort the exclusions */
2591     nra = 0;
2592     for (i = 0; (i < b2->nr); i++)
2593     {
2594         if (b2->nra[i] > 0)
2595         {
2596             /* remove double entries */
2597             qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp);
2598             k = 1;
2599             for (j = 1; (j < b2->nra[i]); j++)
2600             {
2601                 if (b2->a[i][j] != b2->a[i][k-1])
2602                 {
2603                     b2->a[i][k] = b2->a[i][j];
2604                     k++;
2605                 }
2606             }
2607             b2->nra[i] = k;
2608             nra       += b2->nra[i];
2609         }
2610     }
2611     excl->nra = nra;
2612     srenew(excl->a, excl->nra);
2613
2614     b2_to_b(b2, excl);
2615 }
2616
2617 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
2618                            t_nbparam ***nbparam, t_nbparam ***pair)
2619 {
2620     t_atom  atom;
2621     t_param param;
2622     int     i, nr;
2623
2624     /* Define an atom type with all parameters set to zero (no interactions) */
2625     atom.q = 0.0;
2626     atom.m = 0.0;
2627     /* Type for decoupled atoms could be anything,
2628      * this should be changed automatically later when required.
2629      */
2630     atom.ptype = eptAtom;
2631     for (i = 0; (i < MAXFORCEPARAM); i++)
2632     {
2633         param.c[i] = 0.0;
2634     }
2635
2636     nr = add_atomtype(at, symtab, &atom, "decoupled", &param, -1, 0.0, 0.0, 0.0, 0, 0, 0);
2637
2638     /* Add space in the non-bonded parameters matrix */
2639     realloc_nb_params(at, nbparam, pair);
2640
2641     return nr;
2642 }
2643
2644 static void convert_pairs_to_pairsQ(t_params *plist,
2645                                     real fudgeQQ, t_atoms *atoms)
2646 {
2647     t_param *paramp1, *paramp2, *paramnew;
2648     int      i, j, p1nr, p2nr, p2newnr;
2649
2650     /* Add the pair list to the pairQ list */
2651     p1nr    = plist[F_LJ14].nr;
2652     p2nr    = plist[F_LJC14_Q].nr;
2653     p2newnr = p1nr + p2nr;
2654     snew(paramnew, p2newnr);
2655
2656     paramp1             = plist[F_LJ14].param;
2657     paramp2             = plist[F_LJC14_Q].param;
2658
2659     /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2660        it may be possible to just ADD the converted F_LJ14 array
2661        to the old F_LJC14_Q array, but since we have to create
2662        a new sized memory structure, better just to deep copy it all.
2663      */
2664
2665     for (i = 0; i < p2nr; i++)
2666     {
2667         /* Copy over parameters */
2668         for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */
2669         {
2670             paramnew[i].c[j] = paramp2[i].c[j];
2671         }
2672
2673         /* copy over atoms */
2674         for (j = 0; j < 2; j++)
2675         {
2676             paramnew[i].a[j] = paramp2[i].a[j];
2677         }
2678     }
2679
2680     for (i = p2nr; i < p2newnr; i++)
2681     {
2682         j             = i-p2nr;
2683
2684         /* Copy over parameters */
2685         paramnew[i].c[0] = fudgeQQ;
2686         paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q;
2687         paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q;
2688         paramnew[i].c[3] = paramp1[j].c[0];
2689         paramnew[i].c[4] = paramp1[j].c[1];
2690
2691         /* copy over atoms */
2692         paramnew[i].a[0] = paramp1[j].a[0];
2693         paramnew[i].a[1] = paramp1[j].a[1];
2694     }
2695
2696     /* free the old pairlists */
2697     sfree(plist[F_LJC14_Q].param);
2698     sfree(plist[F_LJ14].param);
2699
2700     /* now assign the new data to the F_LJC14_Q structure */
2701     plist[F_LJC14_Q].param   = paramnew;
2702     plist[F_LJC14_Q].nr      = p2newnr;
2703
2704     /* Empty the LJ14 pairlist */
2705     plist[F_LJ14].nr    = 0;
2706     plist[F_LJ14].param = nullptr;
2707 }
2708
2709 static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp, warninp_t wi)
2710 {
2711     int       n, ntype, i, j, k;
2712     t_atom   *atom;
2713     t_blocka *excl;
2714     gmx_bool  bExcl;
2715     t_param   param;
2716     char      errbuf[STRLEN];
2717
2718     n    = mol->atoms.nr;
2719     atom = mol->atoms.atom;
2720
2721     ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->nr)));
2722     GMX_ASSERT(ntype * ntype == nbp->nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
2723
2724     for (i = 0; i < MAXATOMLIST; i++)
2725     {
2726         param.a[i] = NOTSET;
2727     }
2728     for (i = 0; i < MAXFORCEPARAM; i++)
2729     {
2730         param.c[i] = NOTSET;
2731     }
2732
2733     /* Add a pair interaction for all non-excluded atom pairs */
2734     excl = &mol->excls;
2735     for (i = 0; i < n; i++)
2736     {
2737         for (j = i+1; j < n; j++)
2738         {
2739             bExcl = FALSE;
2740             for (k = excl->index[i]; k < excl->index[i+1]; k++)
2741             {
2742                 if (excl->a[k] == j)
2743                 {
2744                     bExcl = TRUE;
2745                 }
2746             }
2747             if (!bExcl)
2748             {
2749                 if (nb_funct != F_LJ)
2750                 {
2751                     sprintf(errbuf, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2752                     warning_error_and_exit(wi, errbuf, FARGS);
2753                 }
2754                 param.a[0] = i;
2755                 param.a[1] = j;
2756                 param.c[0] = atom[i].q;
2757                 param.c[1] = atom[j].q;
2758                 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2759                 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2760                 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], &param);
2761             }
2762         }
2763     }
2764 }
2765
2766 static void set_excl_all(t_blocka *excl)
2767 {
2768     int nat, i, j, k;
2769
2770     /* Get rid of the current exclusions and exclude all atom pairs */
2771     nat       = excl->nr;
2772     excl->nra = nat*nat;
2773     srenew(excl->a, excl->nra);
2774     k = 0;
2775     for (i = 0; i < nat; i++)
2776     {
2777         excl->index[i] = k;
2778         for (j = 0; j < nat; j++)
2779         {
2780             excl->a[k++] = j;
2781         }
2782     }
2783     excl->index[nat] = k;
2784 }
2785
2786 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2787                            int couple_lam0, int couple_lam1,
2788                            const char *mol_name, warninp_t wi)
2789 {
2790     int  i;
2791     char errbuf[STRLEN];
2792
2793     for (i = 0; i < atoms->nr; i++)
2794     {
2795         t_atom *atom;
2796
2797         atom = &atoms->atom[i];
2798
2799         if (atom->qB != atom->q || atom->typeB != atom->type)
2800         {
2801             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.",
2802                     i + 1, mol_name, "couple-moltype");
2803             warning_error_and_exit(wi, errbuf, FARGS);
2804         }
2805
2806         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2807         {
2808             atom->q     = 0.0;
2809         }
2810         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2811         {
2812             atom->type  = atomtype_decouple;
2813         }
2814         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2815         {
2816             atom->qB    = 0.0;
2817         }
2818         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2819         {
2820             atom->typeB = atomtype_decouple;
2821         }
2822     }
2823 }
2824
2825 void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
2826                             int couple_lam0, int couple_lam1,
2827                             gmx_bool bCoupleIntra, int nb_funct, t_params *nbp,
2828                             warninp_t wi)
2829 {
2830     convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2831     if (!bCoupleIntra)
2832     {
2833         generate_LJCpairsNB(mol, nb_funct, nbp, wi);
2834         set_excl_all(&mol->excls);
2835     }
2836     decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1,
2837                    *mol->name, wi);
2838 }