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