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 /* Returns the number of exact atom type matches, i.e. non wild-card matches,
1580  * returns -1 when there are no matches at all.
1581  */
1582 static int natom_match(t_param *pi,
1583                        int type_i, int type_j, int type_k, int type_l,
1584                        const gpp_atomtype_t atype)
1585 {
1586     if ((pi->ai() == -1 || get_atomtype_batype(type_i, atype) == pi->ai()) &&
1587         (pi->aj() == -1 || get_atomtype_batype(type_j, atype) == pi->aj()) &&
1588         (pi->ak() == -1 || get_atomtype_batype(type_k, atype) == pi->ak()) &&
1589         (pi->al() == -1 || get_atomtype_batype(type_l, atype) == pi->al()))
1590     {
1591         return
1592             (pi->ai() == -1 ? 0 : 1) +
1593             (pi->aj() == -1 ? 0 : 1) +
1594             (pi->ak() == -1 ? 0 : 1) +
1595             (pi->al() == -1 ? 0 : 1);
1596     }
1597     else
1598     {
1599         return -1;
1600     }
1601 }
1602
1603 static gmx_bool default_params(int ftype, t_params bt[],
1604                                t_atoms *at, gpp_atomtype_t atype,
1605                                t_param *p, gmx_bool bB,
1606                                t_param **param_def,
1607                                int *nparam_def)
1608 {
1609     int          nparam_found;
1610     gmx_bool     bFound, bSame;
1611     t_param     *pi    = NULL;
1612     t_param     *pj    = NULL;
1613     int          nr    = bt[ftype].nr;
1614     int          nral  = NRAL(ftype);
1615     int          nrfpA = interaction_function[ftype].nrfpA;
1616     int          nrfpB = interaction_function[ftype].nrfpB;
1617
1618     if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1619     {
1620         return TRUE;
1621     }
1622
1623
1624     bFound       = FALSE;
1625     nparam_found = 0;
1626     if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1627     {
1628         int nmatch_max = -1;
1629         int i          = -1;
1630         int t;
1631
1632         /* For dihedrals we allow wildcards. We choose the first type
1633          * that has the most real matches, i.e. non-wildcard matches.
1634          */
1635         for (t = 0; ((t < nr) && nmatch_max < 4); t++)
1636         {
1637             int      nmatch;
1638             t_param *pt;
1639
1640             pt = &(bt[ftype].param[t]);
1641             if (bB)
1642             {
1643                 nmatch = natom_match(pt, at->atom[p->ai()].typeB, at->atom[p->aj()].typeB, at->atom[p->ak()].typeB, at->atom[p->al()].typeB, atype);
1644             }
1645             else
1646             {
1647                 nmatch = natom_match(pt, at->atom[p->ai()].type, at->atom[p->aj()].type, at->atom[p->ak()].type, at->atom[p->al()].type, atype);
1648             }
1649             if (nmatch > nmatch_max)
1650             {
1651                 nmatch_max = nmatch;
1652                 i          = t;
1653                 bFound     = TRUE;
1654             }
1655         }
1656
1657         if (bFound == TRUE)
1658         {
1659             int j;
1660
1661             pi    = &(bt[ftype].param[i]);
1662             nparam_found++;
1663
1664             /* Find additional matches for this dihedral - necessary
1665              * for ftype==9.
1666              * The rule in that case is that additional matches
1667              * HAVE to be on adjacent lines!
1668              */
1669             bSame = TRUE;
1670             /* Continue from current i value */
1671             for (j = i + 2; j < nr && bSame; j += 2)
1672             {
1673                 pj    = &(bt[ftype].param[j]);
1674                 bSame = (pi->ai() == pj->ai() && pi->aj() == pj->aj() && pi->ak() == pj->ak() && pi->al() == pj->al());
1675                 if (bSame)
1676                 {
1677                     nparam_found++;
1678                 }
1679                 /* nparam_found will be increased as long as the numbers match */
1680             }
1681         }
1682     }
1683     else   /* Not a dihedral */
1684     {
1685         int i, j;
1686
1687         for (i = 0; ((i < nr) && !bFound); i++)
1688         {
1689             pi = &(bt[ftype].param[i]);
1690             if (bB)
1691             {
1692                 for (j = 0; ((j < nral) &&
1693                              (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++)
1694                 {
1695                     ;
1696                 }
1697             }
1698             else
1699             {
1700                 for (j = 0; ((j < nral) &&
1701                              (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++)
1702                 {
1703                     ;
1704                 }
1705             }
1706             bFound = (j == nral);
1707         }
1708         if (bFound)
1709         {
1710             nparam_found = 1;
1711         }
1712     }
1713
1714     *param_def  = pi;
1715     *nparam_def = nparam_found;
1716
1717     return bFound;
1718 }
1719
1720
1721
1722 void push_bond(directive d, t_params bondtype[], t_params bond[],
1723                t_atoms *at, gpp_atomtype_t atype, char *line,
1724                gmx_bool bBonded, gmx_bool bGenPairs, real fudgeQQ,
1725                gmx_bool bZero, gmx_bool *bWarn_copy_A_B,
1726                warninp_t wi)
1727 {
1728     const char  *aaformat[MAXATOMLIST] = {
1729         "%d%d",
1730         "%d%d%d",
1731         "%d%d%d%d",
1732         "%d%d%d%d%d",
1733         "%d%d%d%d%d%d",
1734         "%d%d%d%d%d%d%d"
1735     };
1736     const char  *asformat[MAXATOMLIST] = {
1737         "%*s%*s",
1738         "%*s%*s%*s",
1739         "%*s%*s%*s%*s",
1740         "%*s%*s%*s%*s%*s",
1741         "%*s%*s%*s%*s%*s%*s",
1742         "%*s%*s%*s%*s%*s%*s%*s"
1743     };
1744     const char  *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1745     int          nr, i, j, nral, nral_fmt, nread, ftype;
1746     char         format[STRLEN];
1747     /* One force parameter more, so we can check if we read too many */
1748     double       cc[MAXFORCEPARAM+1];
1749     int          aa[MAXATOMLIST+1];
1750     t_param      param, *param_defA, *param_defB;
1751     gmx_bool     bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE;
1752     int          nparam_defA, nparam_defB;
1753     char         errbuf[256];
1754
1755     nparam_defA = nparam_defB = 0;
1756
1757     ftype = ifunc_index(d, 1);
1758     nral  = NRAL(ftype);
1759     for (j = 0; j < MAXATOMLIST; j++)
1760     {
1761         aa[j] = NOTSET;
1762     }
1763     bDef = (NRFP(ftype) > 0);
1764
1765     if (ftype == F_SETTLE)
1766     {
1767         /* SETTLE acts on 3 atoms, but the topology format only specifies
1768          * the first atom (for historical reasons).
1769          */
1770         nral_fmt = 1;
1771     }
1772     else
1773     {
1774         nral_fmt = nral;
1775     }
1776
1777     nread = sscanf(line, aaformat[nral_fmt-1],
1778                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1779
1780     if (ftype == F_SETTLE)
1781     {
1782         aa[3] = aa[1];
1783         aa[1] = aa[0] + 1;
1784         aa[2] = aa[0] + 2;
1785     }
1786
1787     if (nread < nral_fmt)
1788     {
1789         too_few(wi);
1790         return;
1791     }
1792     else if (nread > nral_fmt)
1793     {
1794         /* this is a hack to allow for virtual sites with swapped parity */
1795         bSwapParity = (aa[nral] < 0);
1796         if (bSwapParity)
1797         {
1798             aa[nral] = -aa[nral];
1799         }
1800         ftype = ifunc_index(d, aa[nral]);
1801         if (bSwapParity)
1802         {
1803             switch (ftype)
1804             {
1805                 case F_VSITE3FAD:
1806                 case F_VSITE3OUT:
1807                     break;
1808                 default:
1809                     gmx_fatal(FARGS, "Negative function types only allowed for %s and %s",
1810                               interaction_function[F_VSITE3FAD].longname,
1811                               interaction_function[F_VSITE3OUT].longname);
1812             }
1813         }
1814     }
1815
1816
1817     /* Check for double atoms and atoms out of bounds */
1818     for (i = 0; (i < nral); i++)
1819     {
1820         if (aa[i] < 1 || aa[i] > at->nr)
1821         {
1822             gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
1823                       "Atom index (%d) in %s out of bounds (1-%d).\n"
1824                       "This probably means that you have inserted topology section \"%s\"\n"
1825                       "in a part belonging to a different molecule than you intended to.\n"
1826                       "In that case move the \"%s\" section to the right molecule.",
1827                       get_warning_file(wi), get_warning_line(wi),
1828                       aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
1829         }
1830         for (j = i+1; (j < nral); j++)
1831         {
1832             if (aa[i] == aa[j])
1833             {
1834                 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1835                 warning(wi, errbuf);
1836             }
1837         }
1838     }
1839
1840     /* default force parameters  */
1841     for (j = 0; (j < MAXATOMLIST); j++)
1842     {
1843         param.a[j] = aa[j]-1;
1844     }
1845     for (j = 0; (j < MAXFORCEPARAM); j++)
1846     {
1847         param.c[j] = 0.0;
1848     }
1849
1850     /* Get force params for normal and free energy perturbation
1851      * studies, as determined by types!
1852      */
1853
1854     if (bBonded)
1855     {
1856         bFoundA = default_params(ftype, bondtype, at, atype, &param, FALSE, &param_defA, &nparam_defA);
1857         if (bFoundA)
1858         {
1859             /* Copy the A-state and B-state default parameters. */
1860             GMX_ASSERT(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM, "Bonded interactions may have at most 12 parameters");
1861             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
1862             {
1863                 param.c[j] = param_defA->c[j];
1864             }
1865         }
1866         bFoundB = default_params(ftype, bondtype, at, atype, &param, TRUE, &param_defB, &nparam_defB);
1867         if (bFoundB)
1868         {
1869             /* Copy only the B-state default parameters */
1870             for (j = NRFPA(ftype); (j < NRFP(ftype)); j++)
1871             {
1872                 param.c[j] = param_defB->c[j];
1873             }
1874         }
1875     }
1876     else if (ftype == F_LJ14)
1877     {
1878         bFoundA = default_nb_params(ftype, bondtype, at, &param, 0, FALSE, bGenPairs);
1879         bFoundB = default_nb_params(ftype, bondtype, at, &param, 0, TRUE, bGenPairs);
1880     }
1881     else if (ftype == F_LJC14_Q)
1882     {
1883         param.c[0] = fudgeQQ;
1884         /* Fill in the A-state charges as default parameters */
1885         param.c[1] = at->atom[param.a[0]].q;
1886         param.c[2] = at->atom[param.a[1]].q;
1887         /* The default LJ parameters are the standard 1-4 parameters */
1888         bFoundA = default_nb_params(F_LJ14, bondtype, at, &param, 3, FALSE, bGenPairs);
1889         bFoundB = TRUE;
1890     }
1891     else if (ftype == F_LJC_PAIRS_NB)
1892     {
1893         /* Defaults are not supported here */
1894         bFoundA = FALSE;
1895         bFoundB = TRUE;
1896     }
1897     else
1898     {
1899         gmx_incons("Unknown function type in push_bond");
1900     }
1901
1902     if (nread > nral_fmt)
1903     {
1904         /* Manually specified parameters - in this case we discard multiple torsion info! */
1905
1906         strcpy(format, asformat[nral_fmt-1]);
1907         strcat(format, ccformat);
1908
1909         nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5],
1910                        &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
1911
1912         if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1913         {
1914             /* We only have to issue a warning if these atoms are perturbed! */
1915             bPert = FALSE;
1916             for (j = 0; (j < nral); j++)
1917             {
1918                 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1919             }
1920
1921             if (bPert && *bWarn_copy_A_B)
1922             {
1923                 sprintf(errbuf,
1924                         "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1925                 warning(wi, errbuf);
1926                 *bWarn_copy_A_B = FALSE;
1927             }
1928
1929             /* If only the A parameters were specified, copy them to the B state */
1930             /* The B-state parameters correspond to the first nrfpB
1931              * A-state parameters.
1932              */
1933             for (j = 0; (j < NRFPB(ftype)); j++)
1934             {
1935                 cc[nread++] = cc[j];
1936             }
1937         }
1938
1939         /* If nread was 0 or EOF, no parameters were read => use defaults.
1940          * If nread was nrfpA we copied above so nread=nrfp.
1941          * If nread was nrfp we are cool.
1942          * For F_LJC14_Q we allow supplying fudgeQQ only.
1943          * Anything else is an error!
1944          */
1945         if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1946             !(ftype == F_LJC14_Q && nread == 1))
1947         {
1948             gmx_fatal(FARGS, "Incorrect number of parameters - found %d, expected %d or %d for %s.",
1949                       nread, NRFPA(ftype), NRFP(ftype),
1950                       interaction_function[ftype].longname);
1951         }
1952
1953         for (j = 0; (j < nread); j++)
1954         {
1955             param.c[j] = cc[j];
1956         }
1957
1958         /* Check whether we have to use the defaults */
1959         if (nread == NRFP(ftype))
1960         {
1961             bDef = FALSE;
1962         }
1963     }
1964     else
1965     {
1966         nread = 0;
1967     }
1968     /* nread now holds the number of force parameters read! */
1969
1970     if (bDef)
1971     {
1972         /* Use defaults */
1973         /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1974         if (ftype == F_PDIHS)
1975         {
1976             if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB)))
1977             {
1978                 sprintf(errbuf,
1979                         "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1980                         "Please specify perturbed parameters manually for this torsion in your topology!");
1981                 warning_error(wi, errbuf);
1982             }
1983         }
1984
1985         if (nread > 0 && nread < NRFPA(ftype))
1986         {
1987             /* Issue an error, do not use defaults */
1988             sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
1989             warning_error(wi, errbuf);
1990         }
1991
1992         if (nread == 0 || nread == EOF)
1993         {
1994             if (!bFoundA)
1995             {
1996                 if (interaction_function[ftype].flags & IF_VSITE)
1997                 {
1998                     /* set them to NOTSET, will be calculated later */
1999                     for (j = 0; (j < MAXFORCEPARAM); j++)
2000                     {
2001                         param.c[j] = NOTSET;
2002                     }
2003
2004                     if (bSwapParity)
2005                     {
2006                         param.c1() = -1; /* flag to swap parity of vsite construction */
2007                     }
2008                 }
2009                 else
2010                 {
2011                     if (bZero)
2012                     {
2013                         fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
2014                                 interaction_function[ftype].longname);
2015                     }
2016                     else
2017                     {
2018                         sprintf(errbuf, "No default %s types", interaction_function[ftype].longname);
2019                         warning_error(wi, errbuf);
2020                     }
2021                 }
2022             }
2023             else
2024             {
2025                 if (bSwapParity)
2026                 {
2027                     switch (ftype)
2028                     {
2029                         case F_VSITE3FAD:
2030                             param.c0() = 360 - param.c0();
2031                             break;
2032                         case F_VSITE3OUT:
2033                             param.c2() = -param.c2();
2034                             break;
2035                     }
2036                 }
2037             }
2038             if (!bFoundB)
2039             {
2040                 /* We only have to issue a warning if these atoms are perturbed! */
2041                 bPert = FALSE;
2042                 for (j = 0; (j < nral); j++)
2043                 {
2044                     bPert = bPert || PERTURBED(at->atom[param.a[j]]);
2045                 }
2046
2047                 if (bPert)
2048                 {
2049                     sprintf(errbuf, "No default %s types for perturbed atoms, "
2050                             "using normal values", interaction_function[ftype].longname);
2051                     warning(wi, errbuf);
2052                 }
2053             }
2054         }
2055     }
2056
2057     if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
2058         && param.c[5] != param.c[2])
2059     {
2060         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2061                   "             %s multiplicity can not be perturbed %f!=%f",
2062                   get_warning_file(wi), get_warning_line(wi),
2063                   interaction_function[ftype].longname,
2064                   param.c[2], param.c[5]);
2065     }
2066
2067     if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
2068     {
2069         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2070                   "             %s table number can not be perturbed %d!=%d",
2071                   get_warning_file(wi), get_warning_line(wi),
2072                   interaction_function[ftype].longname,
2073                   (int)(param.c[0]+0.5), (int)(param.c[2]+0.5));
2074     }
2075
2076     /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2077     if (ftype == F_RBDIHS)
2078     {
2079         nr = 0;
2080         for (i = 0; i < NRFP(ftype); i++)
2081         {
2082             if (param.c[i] != 0)
2083             {
2084                 nr++;
2085             }
2086         }
2087         if (nr == 0)
2088         {
2089             return;
2090         }
2091     }
2092
2093     /* Put the values in the appropriate arrays */
2094     add_param_to_list (&bond[ftype], &param);
2095
2096     /* Push additional torsions from FF for ftype==9 if we have them.
2097      * We have already checked that the A/B states do not differ in this case,
2098      * so we do not have to double-check that again, or the vsite stuff.
2099      * In addition, those torsions cannot be automatically perturbed.
2100      */
2101     if (bDef && ftype == F_PDIHS)
2102     {
2103         for (i = 1; i < nparam_defA; i++)
2104         {
2105             /* Advance pointer! */
2106             param_defA += 2;
2107             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
2108             {
2109                 param.c[j] = param_defA->c[j];
2110             }
2111             /* And push the next term for this torsion */
2112             add_param_to_list (&bond[ftype], &param);
2113         }
2114     }
2115 }
2116
2117 void push_cmap(directive d, t_params bondtype[], t_params bond[],
2118                t_atoms *at, gpp_atomtype_t atype, char *line,
2119                warninp_t wi)
2120 {
2121     const char *aaformat[MAXATOMLIST+1] =
2122     {
2123         "%d",
2124         "%d%d",
2125         "%d%d%d",
2126         "%d%d%d%d",
2127         "%d%d%d%d%d",
2128         "%d%d%d%d%d%d",
2129         "%d%d%d%d%d%d%d"
2130     };
2131
2132     int      i, j, ftype, nral, nread, ncmap_params;
2133     int      cmap_type;
2134     int      aa[MAXATOMLIST+1];
2135     char     errbuf[256];
2136     gmx_bool bFound;
2137     t_param  param;
2138
2139     ftype        = ifunc_index(d, 1);
2140     nral         = NRAL(ftype);
2141     ncmap_params = 0;
2142
2143     nread = sscanf(line, aaformat[nral-1],
2144                    &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2145
2146     if (nread < nral)
2147     {
2148         too_few(wi);
2149         return;
2150     }
2151     else if (nread == nral)
2152     {
2153         ftype = ifunc_index(d, 1);
2154     }
2155
2156     /* Check for double atoms and atoms out of bounds */
2157     for (i = 0; i < nral; i++)
2158     {
2159         if (aa[i] < 1 || aa[i] > at->nr)
2160         {
2161             gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2162                       "Atom index (%d) in %s out of bounds (1-%d).\n"
2163                       "This probably means that you have inserted topology section \"%s\"\n"
2164                       "in a part belonging to a different molecule than you intended to.\n"
2165                       "In that case move the \"%s\" section to the right molecule.",
2166                       get_warning_file(wi), get_warning_line(wi),
2167                       aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2168         }
2169
2170         for (j = i+1; (j < nral); j++)
2171         {
2172             if (aa[i] == aa[j])
2173             {
2174                 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2175                 warning(wi, errbuf);
2176             }
2177         }
2178     }
2179
2180     /* default force parameters  */
2181     for (j = 0; (j < MAXATOMLIST); j++)
2182     {
2183         param.a[j] = aa[j]-1;
2184     }
2185     for (j = 0; (j < MAXFORCEPARAM); j++)
2186     {
2187         param.c[j] = 0.0;
2188     }
2189
2190     /* Get the cmap type for this cmap angle */
2191     bFound = default_cmap_params(bondtype, at, atype, &param, FALSE, &cmap_type, &ncmap_params);
2192
2193     /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2194     if (bFound && ncmap_params == 1)
2195     {
2196         /* Put the values in the appropriate arrays */
2197         param.c[0] = cmap_type;
2198         add_param_to_list(&bond[ftype], &param);
2199     }
2200     else
2201     {
2202         /* This is essentially the same check as in default_cmap_params() done one more time */
2203         gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2204                   param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2205     }
2206 }
2207
2208
2209
2210 void push_vsitesn(directive d, t_params bond[],
2211                   t_atoms *at, char *line,
2212                   warninp_t wi)
2213 {
2214     char   *ptr;
2215     int     type, ftype, j, n, ret, nj, a;
2216     int    *atc    = NULL;
2217     double *weight = NULL, weight_tot;
2218     t_param param;
2219
2220     /* default force parameters  */
2221     for (j = 0; (j < MAXATOMLIST); j++)
2222     {
2223         param.a[j] = NOTSET;
2224     }
2225     for (j = 0; (j < MAXFORCEPARAM); j++)
2226     {
2227         param.c[j] = 0.0;
2228     }
2229
2230     ptr  = line;
2231     ret  = sscanf(ptr, "%d%n", &a, &n);
2232     ptr += n;
2233     if (ret == 0)
2234     {
2235         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2236                   "             Expected an atom index in section \"%s\"",
2237                   get_warning_file(wi), get_warning_line(wi),
2238                   dir2str(d));
2239     }
2240
2241     param.a[0] = a - 1;
2242
2243     sscanf(ptr, "%d%n", &type, &n);
2244     ptr  += n;
2245     ftype = ifunc_index(d, type);
2246
2247     weight_tot = 0;
2248     nj         = 0;
2249     do
2250     {
2251         ret  = sscanf(ptr, "%d%n", &a, &n);
2252         ptr += n;
2253         if (ret > 0)
2254         {
2255             if (nj % 20 == 0)
2256             {
2257                 srenew(atc, nj+20);
2258                 srenew(weight, nj+20);
2259             }
2260             atc[nj] = a - 1;
2261             switch (type)
2262             {
2263                 case 1:
2264                     weight[nj] = 1;
2265                     break;
2266                 case 2:
2267                     /* Here we use the A-state mass as a parameter.
2268                      * Note that the B-state mass has no influence.
2269                      */
2270                     weight[nj] = at->atom[atc[nj]].m;
2271                     break;
2272                 case 3:
2273                     weight[nj] = -1;
2274                     ret        = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2275                     ptr       += n;
2276                     if (weight[nj] < 0)
2277                     {
2278                         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2279                                   "             No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
2280                                   get_warning_file(wi), get_warning_line(wi),
2281                                   nj+1, atc[nj]+1);
2282                     }
2283                     break;
2284                 default:
2285                     gmx_fatal(FARGS, "Unknown vsiten type %d", type);
2286             }
2287             weight_tot += weight[nj];
2288             nj++;
2289         }
2290     }
2291     while (ret > 0);
2292
2293     if (nj == 0)
2294     {
2295         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2296                   "             Expected more than one atom index in section \"%s\"",
2297                   get_warning_file(wi), get_warning_line(wi),
2298                   dir2str(d));
2299     }
2300
2301     if (weight_tot == 0)
2302     {
2303         gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2304                   "             The total mass of the construting atoms is zero",
2305                   get_warning_file(wi), get_warning_line(wi));
2306     }
2307
2308     for (j = 0; j < nj; j++)
2309     {
2310         param.a[1] = atc[j];
2311         param.c[0] = nj;
2312         param.c[1] = weight[j]/weight_tot;
2313         /* Put the values in the appropriate arrays */
2314         add_param_to_list (&bond[ftype], &param);
2315     }
2316
2317     sfree(atc);
2318     sfree(weight);
2319 }
2320
2321 void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
2322               int *nrcopies,
2323               warninp_t wi)
2324 {
2325     char type[STRLEN];
2326
2327     if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2328     {
2329         too_few(wi);
2330         return;
2331     }
2332
2333     /* Search moleculename.
2334      * Here we originally only did case insensitive matching. But because
2335      * some PDB files can have many chains and use case to generate more
2336      * chain-identifiers, which in turn end up in our moleculetype name,
2337      * we added support for case-sensitivity.
2338      */
2339     int nrcs    = 0;
2340     int nrci    = 0;
2341     int matchci = -1;
2342     int matchcs = -1;
2343     for (int i = 0; i < nrmols; i++)
2344     {
2345         if (strcmp(type, *(mols[i].name)) == 0)
2346         {
2347             nrcs++;
2348             matchcs = i;
2349         }
2350         if (gmx_strcasecmp(type, *(mols[i].name)) == 0)
2351         {
2352             nrci++;
2353             matchci = i;
2354         }
2355     }
2356
2357     if (nrcs == 1)
2358     {
2359         // select the case sensitive match
2360         *whichmol = matchcs;
2361     }
2362     else
2363     {
2364         // avoid matching case-insensitive when we have multiple matches
2365         if (nrci > 1)
2366         {
2367             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);
2368         }
2369         if (nrci == 1)
2370         {
2371             // select the unique case insensitive match
2372             *whichmol = matchci;
2373         }
2374         else
2375         {
2376             gmx_fatal(FARGS, "No such moleculetype %s", type);
2377         }
2378     }
2379 }
2380
2381 void init_block2(t_block2 *b2, int natom)
2382 {
2383     int i;
2384
2385     b2->nr = natom;
2386     snew(b2->nra, b2->nr);
2387     snew(b2->a, b2->nr);
2388     for (i = 0; (i < b2->nr); i++)
2389     {
2390         b2->a[i] = NULL;
2391     }
2392 }
2393
2394 void done_block2(t_block2 *b2)
2395 {
2396     int i;
2397
2398     if (b2->nr)
2399     {
2400         for (i = 0; (i < b2->nr); i++)
2401         {
2402             sfree(b2->a[i]);
2403         }
2404         sfree(b2->a);
2405         sfree(b2->nra);
2406         b2->nr = 0;
2407     }
2408 }
2409
2410 void push_excl(char *line, t_block2 *b2)
2411 {
2412     int  i, j;
2413     int  n;
2414     char base[STRLEN], format[STRLEN];
2415
2416     if (sscanf(line, "%d", &i) == 0)
2417     {
2418         return;
2419     }
2420
2421     if ((1 <= i) && (i <= b2->nr))
2422     {
2423         i--;
2424     }
2425     else
2426     {
2427         if (debug)
2428         {
2429             fprintf(debug, "Unbound atom %d\n", i-1);
2430         }
2431         return;
2432     }
2433     strcpy(base, "%*d");
2434     do
2435     {
2436         strcpy(format, base);
2437         strcat(format, "%d");
2438         n = sscanf(line, format, &j);
2439         if (n == 1)
2440         {
2441             if ((1 <= j) && (j <= b2->nr))
2442             {
2443                 j--;
2444                 srenew(b2->a[i], ++(b2->nra[i]));
2445                 b2->a[i][b2->nra[i]-1] = j;
2446                 /* also add the reverse exclusion! */
2447                 srenew(b2->a[j], ++(b2->nra[j]));
2448                 b2->a[j][b2->nra[j]-1] = i;
2449                 strcat(base, "%*d");
2450             }
2451             else
2452             {
2453                 gmx_fatal(FARGS, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
2454             }
2455         }
2456     }
2457     while (n == 1);
2458 }
2459
2460 void b_to_b2(t_blocka *b, t_block2 *b2)
2461 {
2462     int     i;
2463     int     j, a;
2464
2465     for (i = 0; (i < b->nr); i++)
2466     {
2467         for (j = b->index[i]; (j < b->index[i+1]); j++)
2468         {
2469             a = b->a[j];
2470             srenew(b2->a[i], ++b2->nra[i]);
2471             b2->a[i][b2->nra[i]-1] = a;
2472         }
2473     }
2474 }
2475
2476 void b2_to_b(t_block2 *b2, t_blocka *b)
2477 {
2478     int     i, nra;
2479     int     j;
2480
2481     nra = 0;
2482     for (i = 0; (i < b2->nr); i++)
2483     {
2484         b->index[i] = nra;
2485         for (j = 0; (j < b2->nra[i]); j++)
2486         {
2487             b->a[nra+j] = b2->a[i][j];
2488         }
2489         nra += b2->nra[i];
2490     }
2491     /* terminate list */
2492     b->index[i] = nra;
2493 }
2494
2495 static int icomp(const void *v1, const void *v2)
2496 {
2497     return (*((int *) v1))-(*((int *) v2));
2498 }
2499
2500 void merge_excl(t_blocka *excl, t_block2 *b2)
2501 {
2502     int     i, k;
2503     int     j;
2504     int     nra;
2505
2506     if (!b2->nr)
2507     {
2508         return;
2509     }
2510     else if (b2->nr != excl->nr)
2511     {
2512         gmx_fatal(FARGS, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2513                   b2->nr, excl->nr);
2514     }
2515     else if (debug)
2516     {
2517         fprintf(debug, "Entering merge_excl\n");
2518     }
2519
2520     /* First copy all entries from excl to b2 */
2521     b_to_b2(excl, b2);
2522
2523     /* Count and sort the exclusions */
2524     nra = 0;
2525     for (i = 0; (i < b2->nr); i++)
2526     {
2527         if (b2->nra[i] > 0)
2528         {
2529             /* remove double entries */
2530             qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp);
2531             k = 1;
2532             for (j = 1; (j < b2->nra[i]); j++)
2533             {
2534                 if (b2->a[i][j] != b2->a[i][k-1])
2535                 {
2536                     b2->a[i][k] = b2->a[i][j];
2537                     k++;
2538                 }
2539             }
2540             b2->nra[i] = k;
2541             nra       += b2->nra[i];
2542         }
2543     }
2544     excl->nra = nra;
2545     srenew(excl->a, excl->nra);
2546
2547     b2_to_b(b2, excl);
2548 }
2549
2550 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
2551                            t_nbparam ***nbparam, t_nbparam ***pair)
2552 {
2553     t_atom  atom;
2554     t_param param;
2555     int     i, nr;
2556
2557     /* Define an atom type with all parameters set to zero (no interactions) */
2558     atom.q = 0.0;
2559     atom.m = 0.0;
2560     /* Type for decoupled atoms could be anything,
2561      * this should be changed automatically later when required.
2562      */
2563     atom.ptype = eptAtom;
2564     for (i = 0; (i < MAXFORCEPARAM); i++)
2565     {
2566         param.c[i] = 0.0;
2567     }
2568
2569     nr = add_atomtype(at, symtab, &atom, "decoupled", &param, -1, 0.0, 0.0, 0.0, 0, 0, 0);
2570
2571     /* Add space in the non-bonded parameters matrix */
2572     realloc_nb_params(at, nbparam, pair);
2573
2574     return nr;
2575 }
2576
2577 static void convert_pairs_to_pairsQ(t_params *plist,
2578                                     real fudgeQQ, t_atoms *atoms)
2579 {
2580     t_param *paramp1, *paramp2, *paramnew;
2581     int      i, j, p1nr, p2nr, p2newnr;
2582
2583     /* Add the pair list to the pairQ list */
2584     p1nr    = plist[F_LJ14].nr;
2585     p2nr    = plist[F_LJC14_Q].nr;
2586     p2newnr = p1nr + p2nr;
2587     snew(paramnew, p2newnr);
2588
2589     paramp1             = plist[F_LJ14].param;
2590     paramp2             = plist[F_LJC14_Q].param;
2591
2592     /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2593        it may be possible to just ADD the converted F_LJ14 array
2594        to the old F_LJC14_Q array, but since we have to create
2595        a new sized memory structure, better just to deep copy it all.
2596      */
2597
2598     for (i = 0; i < p2nr; i++)
2599     {
2600         /* Copy over parameters */
2601         for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */
2602         {
2603             paramnew[i].c[j] = paramp2[i].c[j];
2604         }
2605
2606         /* copy over atoms */
2607         for (j = 0; j < 2; j++)
2608         {
2609             paramnew[i].a[j] = paramp2[i].a[j];
2610         }
2611     }
2612
2613     for (i = p2nr; i < p2newnr; i++)
2614     {
2615         j             = i-p2nr;
2616
2617         /* Copy over parameters */
2618         paramnew[i].c[0] = fudgeQQ;
2619         paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q;
2620         paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q;
2621         paramnew[i].c[3] = paramp1[j].c[0];
2622         paramnew[i].c[4] = paramp1[j].c[1];
2623
2624         /* copy over atoms */
2625         paramnew[i].a[0] = paramp1[j].a[0];
2626         paramnew[i].a[1] = paramp1[j].a[1];
2627     }
2628
2629     /* free the old pairlists */
2630     sfree(plist[F_LJC14_Q].param);
2631     sfree(plist[F_LJ14].param);
2632
2633     /* now assign the new data to the F_LJC14_Q structure */
2634     plist[F_LJC14_Q].param   = paramnew;
2635     plist[F_LJC14_Q].nr      = p2newnr;
2636
2637     /* Empty the LJ14 pairlist */
2638     plist[F_LJ14].nr    = 0;
2639     plist[F_LJ14].param = NULL;
2640 }
2641
2642 static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp)
2643 {
2644     int       n, ntype, i, j, k;
2645     t_atom   *atom;
2646     t_blocka *excl;
2647     gmx_bool  bExcl;
2648     t_param   param;
2649
2650     n    = mol->atoms.nr;
2651     atom = mol->atoms.atom;
2652
2653     ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->nr)));
2654     GMX_ASSERT(ntype * ntype == nbp->nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
2655
2656     for (i = 0; i < MAXATOMLIST; i++)
2657     {
2658         param.a[i] = NOTSET;
2659     }
2660     for (i = 0; i < MAXFORCEPARAM; i++)
2661     {
2662         param.c[i] = NOTSET;
2663     }
2664
2665     /* Add a pair interaction for all non-excluded atom pairs */
2666     excl = &mol->excls;
2667     for (i = 0; i < n; i++)
2668     {
2669         for (j = i+1; j < n; j++)
2670         {
2671             bExcl = FALSE;
2672             for (k = excl->index[i]; k < excl->index[i+1]; k++)
2673             {
2674                 if (excl->a[k] == j)
2675                 {
2676                     bExcl = TRUE;
2677                 }
2678             }
2679             if (!bExcl)
2680             {
2681                 if (nb_funct != F_LJ)
2682                 {
2683                     gmx_fatal(FARGS, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2684                 }
2685                 param.a[0] = i;
2686                 param.a[1] = j;
2687                 param.c[0] = atom[i].q;
2688                 param.c[1] = atom[j].q;
2689                 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2690                 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2691                 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], &param);
2692             }
2693         }
2694     }
2695 }
2696
2697 static void set_excl_all(t_blocka *excl)
2698 {
2699     int nat, i, j, k;
2700
2701     /* Get rid of the current exclusions and exclude all atom pairs */
2702     nat       = excl->nr;
2703     excl->nra = nat*nat;
2704     srenew(excl->a, excl->nra);
2705     k = 0;
2706     for (i = 0; i < nat; i++)
2707     {
2708         excl->index[i] = k;
2709         for (j = 0; j < nat; j++)
2710         {
2711             excl->a[k++] = j;
2712         }
2713     }
2714     excl->index[nat] = k;
2715 }
2716
2717 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2718                            int couple_lam0, int couple_lam1,
2719                            const char *mol_name)
2720 {
2721     int i;
2722
2723     for (i = 0; i < atoms->nr; i++)
2724     {
2725         t_atom *atom;
2726
2727         atom = &atoms->atom[i];
2728
2729         if (atom->qB != atom->q || atom->typeB != atom->type)
2730         {
2731             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.",
2732                       i + 1, mol_name, "couple-moltype");
2733         }
2734
2735         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2736         {
2737             atom->q     = 0.0;
2738         }
2739         if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2740         {
2741             atom->type  = atomtype_decouple;
2742         }
2743         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2744         {
2745             atom->qB    = 0.0;
2746         }
2747         if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2748         {
2749             atom->typeB = atomtype_decouple;
2750         }
2751     }
2752 }
2753
2754 void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
2755                             int couple_lam0, int couple_lam1,
2756                             gmx_bool bCoupleIntra, int nb_funct, t_params *nbp)
2757 {
2758     convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2759     if (!bCoupleIntra)
2760     {
2761         generate_LJCpairsNB(mol, nb_funct, nbp);
2762         set_excl_all(&mol->excls);
2763     }
2764     decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1,
2765                    *mol->name);
2766 }