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