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