src/gromacs/gmxpreprocess/toppush.cpp

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