src/gromacs/gmxpreprocess/toppush.c

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