src/gromacs/gmxpreprocess/pdb2top.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, 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 <stdio.h>
  42 #include <math.h>
  43 #include <ctype.h>
  44
  45 #include "gromacs/math/vec.h"
  46 #include "macros.h"
  47 #include "gromacs/utility/futil.h"
  48 #include "pdb2top.h"
  49 #include "gpp_nextnb.h"
  50 #include "topdirs.h"
  51 #include "toputil.h"
  52 #include "h_db.h"
  53 #include "pgutil.h"
  54 #include "resall.h"
  55 #include "topio.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/fileio/pdbio.h"
  58 #include "gen_ad.h"
  59 #include "gromacs/fileio/filenm.h"
  60 #include "gen_vsite.h"
  61 #include "add_par.h"
  62 #include "toputil.h"
  63 #include "fflibutil.h"
  64 #include "copyrite.h"
  65
  66 #include "gromacs/fileio/strdb.h"
  67 #include "gromacs/topology/residuetypes.h"
  68 #include "gromacs/topology/symtab.h"
  69 #include "gromacs/utility/exceptions.h"
  70 #include "gromacs/utility/fatalerror.h"
  71 #include "gromacs/utility/programcontext.h"
  72 #include "gromacs/utility/smalloc.h"
  73
  74 /* this must correspond to enum in pdb2top.h */
  75 const char *hh[ehisNR]   = { "HISD", "HISE", "HISH", "HIS1" };
  76
  77 static int missing_atoms(t_restp *rp, int resind, t_atoms *at, int i0, int i)
  78 {
  79     int      j, k, nmiss;
  80     char    *name;
  81     gmx_bool bFound;
  82
  83     nmiss = 0;
  84     for (j = 0; j < rp->natom; j++)
  85     {
  86         name   = *(rp->atomname[j]);
  87         bFound = FALSE;
  88         for (k = i0; k < i; k++)
  89         {
  90             bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]), name));
  91         }
  92         if (!bFound)
  93         {
  94             nmiss++;
  95             fprintf(stderr, "\nWARNING: "
  96                     "atom %s is missing in residue %s %d in the pdb file\n",
  97                     name, *(at->resinfo[resind].name), at->resinfo[resind].nr);
  98             if (name[0] == 'H' || name[0] == 'h')
  99             {
 100                 fprintf(stderr, "         You might need to add atom %s to the hydrogen database of building block %s\n"
 101                         "         in the file %s.hdb (see the manual)\n",
 102                         name, *(at->resinfo[resind].rtp), rp->filebase);
 103             }
 104             fprintf(stderr, "\n");
 105         }
 106     }
 107
 108     return nmiss;
 109 }
 110
 111 gmx_bool is_int(double x)
 112 {
 113     const double tol = 1e-4;
 114     int          ix;
 115
 116     if (x < 0)
 117     {
 118         x = -x;
 119     }
 120     ix = gmx_nint(x);
 121
 122     return (fabs(x-ix) < tol);
 123 }
 124
 125 static void swap_strings(char **s, int i, int j)
 126 {
 127     char *tmp;
 128
 129     tmp  = s[i];
 130     s[i] = s[j];
 131     s[j] = tmp;
 132 }
 133
 134 void
 135 choose_ff(const char *ffsel,
 136           char *forcefield, int ff_maxlen,
 137           char *ffdir, int ffdir_maxlen)
 138 {
 139     int    nff;
 140     char **ffdirs, **ffs, **ffs_dir, *ptr;
 141     int    i, j, sel, cwdsel, nfound;
 142     char   buf[STRLEN], **desc;
 143     FILE  *fp;
 144     char  *pret;
 145
 146     nff = fflib_search_file_in_dirend(fflib_forcefield_itp(),
 147                                       fflib_forcefield_dir_ext(),
 148                                       &ffdirs);
 149
 150     if (nff == 0)
 151     {
 152         gmx_fatal(FARGS, "No force fields found (files with name '%s' in subdirectories ending on '%s')",
 153                   fflib_forcefield_itp(), fflib_forcefield_dir_ext());
 154     }
 155
 156     /* Replace with unix path separators */
 157     if (DIR_SEPARATOR != '/')
 158     {
 159         for (i = 0; i < nff; i++)
 160         {
 161             while ( (ptr = strchr(ffdirs[i], DIR_SEPARATOR)) != NULL)
 162             {
 163                 *ptr = '/';
 164             }
 165         }
 166     }
 167
 168     /* Store the force field names in ffs */
 169     snew(ffs, nff);
 170     snew(ffs_dir, nff);
 171     for (i = 0; i < nff; i++)
 172     {
 173         /* Remove the path from the ffdir name - use our unix standard here! */
 174         ptr = strrchr(ffdirs[i], '/');
 175         if (ptr == NULL)
 176         {
 177             ffs[i]     = strdup(ffdirs[i]);
 178             ffs_dir[i] = low_gmxlibfn(ffdirs[i], FALSE, FALSE);
 179             if (ffs_dir[i] == NULL)
 180             {
 181                 gmx_fatal(FARGS, "Can no longer find file '%s'", ffdirs[i]);
 182             }
 183         }
 184         else
 185         {
 186             ffs[i]     = strdup(ptr+1);
 187             ffs_dir[i] = strdup(ffdirs[i]);
 188         }
 189         ffs_dir[i][strlen(ffs_dir[i])-strlen(ffs[i])-1] = '\0';
 190         /* Remove the extension from the ffdir name */
 191         ffs[i][strlen(ffs[i])-strlen(fflib_forcefield_dir_ext())] = '\0';
 192     }
 193
 194     if (ffsel != NULL)
 195     {
 196         sel     = -1;
 197         cwdsel  = -1;
 198         nfound  = 0;
 199         for (i = 0; i < nff; i++)
 200         {
 201             if (strcmp(ffs[i], ffsel) == 0)
 202             {
 203                 /* Matching ff name */
 204                 sel = i;
 205                 nfound++;
 206
 207                 if (strncmp(ffs_dir[i], ".", 1) == 0)
 208                 {
 209                     cwdsel = i;
 210                 }
 211             }
 212         }
 213
 214         if (cwdsel != -1)
 215         {
 216             sel = cwdsel;
 217         }
 218
 219         if (nfound > 1)
 220         {
 221             if (cwdsel != -1)
 222             {
 223                 fprintf(stderr,
 224                         "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
 225                         "current directory. Use interactive selection (not the -ff option) if\n"
 226                         "you would prefer a different one.\n", ffsel, nfound);
 227             }
 228             else
 229             {
 230                 gmx_fatal(FARGS,
 231                           "Force field '%s' occurs in %d places, but not in the current directory.\n"
 232                           "Run without the -ff switch and select the force field interactively.", ffsel, nfound);
 233             }
 234         }
 235         else if (nfound == 0)
 236         {
 237             gmx_fatal(FARGS, "Could not find force field '%s' in current directory, install tree or GMXDATA path.", ffsel);
 238         }
 239     }
 240     else if (nff > 1)
 241     {
 242         snew(desc, nff);
 243         for (i = 0; (i < nff); i++)
 244         {
 245             sprintf(buf, "%s%c%s%s%c%s",
 246                     ffs_dir[i], DIR_SEPARATOR,
 247                     ffs[i], fflib_forcefield_dir_ext(), DIR_SEPARATOR,
 248                     fflib_forcefield_doc());
 249             if (gmx_fexist(buf))
 250             {
 251                 /* We don't use fflib_open, because we don't want printf's */
 252                 fp = gmx_ffopen(buf, "r");
 253                 snew(desc[i], STRLEN);
 254                 get_a_line(fp, desc[i], STRLEN);
 255                 gmx_ffclose(fp);
 256             }
 257             else
 258             {
 259                 desc[i] = strdup(ffs[i]);
 260             }
 261         }
 262         /* Order force fields from the same dir alphabetically
 263          * and put deprecated force fields at the end.
 264          */
 265         for (i = 0; (i < nff); i++)
 266         {
 267             for (j = i+1; (j < nff); j++)
 268             {
 269                 if (strcmp(ffs_dir[i], ffs_dir[j]) == 0 &&
 270                     ((desc[i][0] == '[' && desc[j][0] != '[') ||
 271                      ((desc[i][0] == '[' || desc[j][0] != '[') &&
 272                       gmx_strcasecmp(desc[i], desc[j]) > 0)))
 273                 {
 274                     swap_strings(ffdirs, i, j);
 275                     swap_strings(ffs, i, j);
 276                     swap_strings(desc, i, j);
 277                 }
 278             }
 279         }
 280
 281         printf("\nSelect the Force Field:\n");
 282         for (i = 0; (i < nff); i++)
 283         {
 284             if (i == 0 || strcmp(ffs_dir[i-1], ffs_dir[i]) != 0)
 285             {
 286                 if (strcmp(ffs_dir[i], ".") == 0)
 287                 {
 288                     printf("From current directory:\n");
 289                 }
 290                 else
 291                 {
 292                     printf("From '%s':\n", ffs_dir[i]);
 293                 }
 294             }
 295             printf("%2d: %s\n", i+1, desc[i]);
 296             sfree(desc[i]);
 297         }
 298         sfree(desc);
 299
 300         sel = -1;
 301         do
 302         {
 303             pret = fgets(buf, STRLEN, stdin);
 304
 305             if (pret != NULL)
 306             {
 307                 sel = strtol(buf, NULL, 10);
 308                 sel--;
 309             }
 310         }
 311         while (pret == NULL || (sel < 0) || (sel >= nff));
 312
 313         /* Check for a current limitation of the fflib code.
 314          * It will always read from the first ff directory in the list.
 315          * This check assumes that the order of ffs matches the order
 316          * in which fflib_open searches ff library files.
 317          */
 318         for (i = 0; i < sel; i++)
 319         {
 320             if (strcmp(ffs[i], ffs[sel]) == 0)
 321             {
 322                 gmx_fatal(FARGS, "Can only select the first of multiple force field entries with directory name '%s%s' in the list. If you want to use the next entry, run pdb2gmx in a different directory or rename or move the force field directory present in the current working directory.",
 323                           ffs[sel], fflib_forcefield_dir_ext());
 324             }
 325         }
 326     }
 327     else
 328     {
 329         sel = 0;
 330     }
 331
 332     if (strlen(ffs[sel]) >= (size_t)ff_maxlen)
 333     {
 334         gmx_fatal(FARGS, "Length of force field name (%d) >= maxlen (%d)",
 335                   strlen(ffs[sel]), ff_maxlen);
 336     }
 337     strcpy(forcefield, ffs[sel]);
 338
 339     if (strlen(ffdirs[sel]) >= (size_t)ffdir_maxlen)
 340     {
 341         gmx_fatal(FARGS, "Length of force field dir (%d) >= maxlen (%d)",
 342                   strlen(ffdirs[sel]), ffdir_maxlen);
 343     }
 344     strcpy(ffdir, ffdirs[sel]);
 345
 346     for (i = 0; (i < nff); i++)
 347     {
 348         sfree(ffdirs[i]);
 349         sfree(ffs[i]);
 350         sfree(ffs_dir[i]);
 351     }
 352     sfree(ffdirs);
 353     sfree(ffs);
 354     sfree(ffs_dir);
 355 }
 356
 357 void choose_watermodel(const char *wmsel, const char *ffdir,
 358                        char **watermodel)
 359 {
 360     const char *fn_watermodels = "watermodels.dat";
 361     char        fn_list[STRLEN];
 362     FILE       *fp;
 363     char        buf[STRLEN];
 364     int         nwm, sel, i;
 365     char      **model;
 366     char       *pret;
 367
 368     if (strcmp(wmsel, "none") == 0)
 369     {
 370         *watermodel = NULL;
 371
 372         return;
 373     }
 374     else if (strcmp(wmsel, "select") != 0)
 375     {
 376         *watermodel = strdup(wmsel);
 377
 378         return;
 379     }
 380
 381     sprintf(fn_list, "%s%c%s", ffdir, DIR_SEPARATOR, fn_watermodels);
 382
 383     if (!fflib_fexist(fn_list))
 384     {
 385         fprintf(stderr, "No file '%s' found, will not include a water model\n",
 386                 fn_watermodels);
 387         *watermodel = NULL;
 388
 389         return;
 390     }
 391
 392     fp = fflib_open(fn_list);
 393     printf("\nSelect the Water Model:\n");
 394     nwm   = 0;
 395     model = NULL;
 396     while (get_a_line(fp, buf, STRLEN))
 397     {
 398         srenew(model, nwm+1);
 399         snew(model[nwm], STRLEN);
 400         sscanf(buf, "%s%n", model[nwm], &i);
 401         if (i > 0)
 402         {
 403             ltrim(buf+i);
 404             fprintf(stderr, "%2d: %s\n", nwm+1, buf+i);
 405             nwm++;
 406         }
 407         else
 408         {
 409             sfree(model[nwm]);
 410         }
 411     }
 412     gmx_ffclose(fp);
 413     fprintf(stderr, "%2d: %s\n", nwm+1, "None");
 414
 415     sel = -1;
 416     do
 417     {
 418         pret = fgets(buf, STRLEN, stdin);
 419
 420         if (pret != NULL)
 421         {
 422             sel = strtol(buf, NULL, 10);
 423             sel--;
 424         }
 425     }
 426     while (pret == NULL || sel < 0 || sel > nwm);
 427
 428     if (sel == nwm)
 429     {
 430         *watermodel = NULL;
 431     }
 432     else
 433     {
 434         *watermodel = strdup(model[sel]);
 435     }
 436
 437     for (i = 0; i < nwm; i++)
 438     {
 439         sfree(model[i]);
 440     }
 441     sfree(model);
 442 }
 443
 444 static int name2type(t_atoms *at, int **cgnr,
 445                      t_restp restp[], gmx_residuetype_t *rt)
 446 {
 447     int         i, j, prevresind, resind, i0, prevcg, cg, curcg;
 448     char       *name;
 449     gmx_bool    bNterm;
 450     double      qt;
 451     int         nmissat;
 452
 453     nmissat = 0;
 454
 455     resind = -1;
 456     bNterm = FALSE;
 457     i0     = 0;
 458     snew(*cgnr, at->nr);
 459     qt     = 0;
 460     prevcg = NOTSET;
 461     curcg  = 0;
 462     cg     = -1;
 463     j      = NOTSET;
 464
 465     for (i = 0; (i < at->nr); i++)
 466     {
 467         prevresind = resind;
 468         if (at->atom[i].resind != resind)
 469         {
 470             gmx_bool bProt;
 471             resind = at->atom[i].resind;
 472             bProt  = gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name));
 473             bNterm = bProt && (resind == 0);
 474             if (resind > 0)
 475             {
 476                 nmissat += missing_atoms(&restp[prevresind], prevresind, at, i0, i);
 477             }
 478             i0 = i;
 479         }
 480         if (at->atom[i].m == 0)
 481         {
 482             if (debug)
 483             {
 484                 fprintf(debug, "atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
 485                         i+1, *(at->atomname[i]), curcg, prevcg,
 486                         j == NOTSET ? NOTSET : restp[resind].cgnr[j]);
 487             }
 488             qt               = 0;
 489             prevcg           = cg;
 490             name             = *(at->atomname[i]);
 491             j                = search_jtype(&restp[resind], name, bNterm);
 492             at->atom[i].type = restp[resind].atom[j].type;
 493             at->atom[i].q    = restp[resind].atom[j].q;
 494             at->atom[i].m    = restp[resind].atom[j].m;
 495             cg               = restp[resind].cgnr[j];
 496             /* A charge group number -1 signals a separate charge group
 497              * for this atom.
 498              */
 499             if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) )
 500             {
 501                 curcg++;
 502             }
 503         }
 504         else
 505         {
 506             if (debug)
 507             {
 508                 fprintf(debug, "atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
 509                         i+1, *(at->atomname[i]), curcg, qt, is_int(qt));
 510             }
 511             cg = -1;
 512             if (is_int(qt))
 513             {
 514                 qt = 0;
 515                 curcg++;
 516             }
 517             qt += at->atom[i].q;
 518         }
 519         (*cgnr)[i]        = curcg;
 520         at->atom[i].typeB = at->atom[i].type;
 521         at->atom[i].qB    = at->atom[i].q;
 522         at->atom[i].mB    = at->atom[i].m;
 523     }
 524     nmissat += missing_atoms(&restp[resind], resind, at, i0, i);
 525
 526     return nmissat;
 527 }
 528
 529 static void print_top_heavy_H(FILE *out, real mHmult)
 530 {
 531     if (mHmult == 2.0)
 532     {
 533         fprintf(out, "; Using deuterium instead of hydrogen\n\n");
 534     }
 535     else if (mHmult == 4.0)
 536     {
 537         fprintf(out, "#define HEAVY_H\n\n");
 538     }
 539     else if (mHmult != 1.0)
 540     {
 541         fprintf(stderr, "WARNING: unsupported proton mass multiplier (%g) "
 542                 "in pdb2top\n", mHmult);
 543     }
 544 }
 545
 546 void print_top_comment(FILE       *out,
 547                        const char *filename,
 548                        const char *ffdir,
 549                        gmx_bool    bITP)
 550 {
 551     char  ffdir_parent[STRLEN];
 552     char *p;
 553
 554     nice_header(out, filename);
 555     fprintf(out, ";\tThis is a %s topology file\n;\n", bITP ? "include" : "standalone");
 556     try
 557     {
 558         gmx::BinaryInformationSettings settings;
 559         settings.generatedByHeader(true);
 560         settings.linePrefix(";\t");
 561         gmx::printBinaryInformation(out, gmx::getProgramContext(), settings);
 562     }
 563     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 564
 565     if (strchr(ffdir, '/') == NULL)
 566     {
 567         fprintf(out, ";\tForce field was read from the standard Gromacs share directory.\n;\n\n");
 568     }
 569     else if (ffdir[0] == '.')
 570     {
 571         fprintf(out, ";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
 572     }
 573     else
 574     {
 575         strncpy(ffdir_parent, ffdir, STRLEN-1);
 576         ffdir_parent[STRLEN-1] = '\0'; /*make sure it is 0-terminated even for long string*/
 577         p = strrchr(ffdir_parent, '/');
 578
 579         *p = '\0';
 580
 581         fprintf(out,
 582                 ";\tForce field data was read from:\n"
 583                 ";\t%s\n"
 584                 ";\n"
 585                 ";\tNote:\n"
 586                 ";\tThis might be a non-standard force field location. When you use this topology, the\n"
 587                 ";\tforce field must either be present in the current directory, or the location\n"
 588                 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
 589                 ffdir_parent);
 590     }
 591 }
 592
 593 void print_top_header(FILE *out, const char *filename,
 594                       gmx_bool bITP, const char *ffdir, real mHmult)
 595 {
 596     const char *p;
 597
 598     print_top_comment(out, filename, ffdir, bITP);
 599
 600     print_top_heavy_H(out, mHmult);
 601     fprintf(out, "; Include forcefield parameters\n");
 602
 603     p = strrchr(ffdir, '/');
 604     p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
 605
 606     fprintf(out, "#include \"%s/%s\"\n\n", p, fflib_forcefield_itp());
 607 }
 608
 609 static void print_top_posre(FILE *out, const char *pr)
 610 {
 611     fprintf(out, "; Include Position restraint file\n");
 612     fprintf(out, "#ifdef POSRES\n");
 613     fprintf(out, "#include \"%s\"\n", pr);
 614     fprintf(out, "#endif\n\n");
 615 }
 616
 617 static void print_top_water(FILE *out, const char *ffdir, const char *water)
 618 {
 619     const char *p;
 620     char        buf[STRLEN];
 621
 622     fprintf(out, "; Include water topology\n");
 623
 624     p = strrchr(ffdir, '/');
 625     p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
 626     fprintf(out, "#include \"%s/%s.itp\"\n", p, water);
 627
 628     fprintf(out, "\n");
 629     fprintf(out, "#ifdef POSRES_WATER\n");
 630     fprintf(out, "; Position restraint for each water oxygen\n");
 631     fprintf(out, "[ position_restraints ]\n");
 632     fprintf(out, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
 633     fprintf(out, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
 634     fprintf(out, "#endif\n");
 635     fprintf(out, "\n");
 636
 637     sprintf(buf, "%s/ions.itp", p);
 638
 639     if (fflib_fexist(buf))
 640     {
 641         fprintf(out, "; Include topology for ions\n");
 642         fprintf(out, "#include \"%s\"\n", buf);
 643         fprintf(out, "\n");
 644     }
 645 }
 646
 647 static void print_top_system(FILE *out, const char *title)
 648 {
 649     fprintf(out, "[ %s ]\n", dir2str(d_system));
 650     fprintf(out, "; Name\n");
 651     fprintf(out, "%s\n\n", title[0] ? title : "Protein");
 652 }
 653
 654 void print_top_mols(FILE *out,
 655                     const char *title, const char *ffdir, const char *water,
 656                     int nincl, char **incls, int nmol, t_mols *mols)
 657 {
 658     int   i;
 659     char *incl;
 660
 661     if (nincl > 0)
 662     {
 663         fprintf(out, "; Include chain topologies\n");
 664         for (i = 0; (i < nincl); i++)
 665         {
 666             incl = strrchr(incls[i], DIR_SEPARATOR);
 667             if (incl == NULL)
 668             {
 669                 incl = incls[i];
 670             }
 671             else
 672             {
 673                 /* Remove the path from the include name */
 674                 incl = incl + 1;
 675             }
 676             fprintf(out, "#include \"%s\"\n", incl);
 677         }
 678         fprintf(out, "\n");
 679     }
 680
 681     if (water)
 682     {
 683         print_top_water(out, ffdir, water);
 684     }
 685     print_top_system(out, title);
 686
 687     if (nmol)
 688     {
 689         fprintf(out, "[ %s ]\n", dir2str(d_molecules));
 690         fprintf(out, "; %-15s %5s\n", "Compound", "#mols");
 691         for (i = 0; (i < nmol); i++)
 692         {
 693             fprintf(out, "%-15s %5d\n", mols[i].name, mols[i].nr);
 694         }
 695     }
 696 }
 697
 698 void write_top(FILE *out, char *pr, char *molname,
 699                t_atoms *at, gmx_bool bRTPresname,
 700                int bts[], t_params plist[], t_excls excls[],
 701                gpp_atomtype_t atype, int *cgnr, int nrexcl)
 702 /* NOTE: nrexcl is not the size of *excl! */
 703 {
 704     if (at && atype && cgnr)
 705     {
 706         fprintf(out, "[ %s ]\n", dir2str(d_moleculetype));
 707         fprintf(out, "; %-15s %5s\n", "Name", "nrexcl");
 708         fprintf(out, "%-15s %5d\n\n", molname ? molname : "Protein", nrexcl);
 709
 710         print_atoms(out, atype, at, cgnr, bRTPresname);
 711         print_bondeds(out, at->nr, d_bonds,      F_BONDS,    bts[ebtsBONDS], plist);
 712         print_bondeds(out, at->nr, d_constraints, F_CONSTR,   0,              plist);
 713         print_bondeds(out, at->nr, d_constraints, F_CONSTRNC, 0,              plist);
 714         print_bondeds(out, at->nr, d_pairs,      F_LJ14,     0,              plist);
 715         print_excl(out, at->nr, excls);
 716         print_bondeds(out, at->nr, d_angles,     F_ANGLES,   bts[ebtsANGLES], plist);
 717         print_bondeds(out, at->nr, d_dihedrals,  F_PDIHS,    bts[ebtsPDIHS], plist);
 718         print_bondeds(out, at->nr, d_dihedrals,  F_IDIHS,    bts[ebtsIDIHS], plist);
 719         print_bondeds(out, at->nr, d_cmap,       F_CMAP,     bts[ebtsCMAP],  plist);
 720         print_bondeds(out, at->nr, d_polarization, F_POLARIZATION,   0,       plist);
 721         print_bondeds(out, at->nr, d_thole_polarization, F_THOLE_POL, 0,       plist);
 722         print_bondeds(out, at->nr, d_vsites2,    F_VSITE2,   0,              plist);
 723         print_bondeds(out, at->nr, d_vsites3,    F_VSITE3,   0,              plist);
 724         print_bondeds(out, at->nr, d_vsites3,    F_VSITE3FD, 0,              plist);
 725         print_bondeds(out, at->nr, d_vsites3,    F_VSITE3FAD, 0,              plist);
 726         print_bondeds(out, at->nr, d_vsites3,    F_VSITE3OUT, 0,              plist);
 727         print_bondeds(out, at->nr, d_vsites4,    F_VSITE4FD, 0,              plist);
 728         print_bondeds(out, at->nr, d_vsites4,    F_VSITE4FDN, 0,             plist);
 729
 730         if (pr)
 731         {
 732             print_top_posre(out, pr);
 733         }
 734     }
 735 }
 736
 737 static atom_id search_res_atom(const char *type, int resind,
 738                                t_atoms *atoms,
 739                                const char *bondtype, gmx_bool bAllowMissing)
 740 {
 741     int i;
 742
 743     for (i = 0; (i < atoms->nr); i++)
 744     {
 745         if (atoms->atom[i].resind == resind)
 746         {
 747             return search_atom(type, i, atoms, bondtype, bAllowMissing);
 748         }
 749     }
 750
 751     return NO_ATID;
 752 }
 753
 754 static void do_ssbonds(t_params *ps, t_atoms *atoms,
 755                        int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
 756 {
 757     int     i, ri, rj;
 758     atom_id ai, aj;
 759
 760     for (i = 0; (i < nssbonds); i++)
 761     {
 762         ri = ssbonds[i].res1;
 763         rj = ssbonds[i].res2;
 764         ai = search_res_atom(ssbonds[i].a1, ri, atoms,
 765                              "special bond", bAllowMissing);
 766         aj = search_res_atom(ssbonds[i].a2, rj, atoms,
 767                              "special bond", bAllowMissing);
 768         if ((ai == NO_ATID) || (aj == NO_ATID))
 769         {
 770             gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
 771                       ssbonds[i].a1, ssbonds[i].a2);
 772         }
 773         add_param(ps, ai, aj, NULL, NULL);
 774     }
 775 }
 776
 777 static void at2bonds(t_params *psb, t_hackblock *hb,
 778                      t_atoms *atoms,
 779                      rvec x[],
 780                      real long_bond_dist, real short_bond_dist)
 781 {
 782     int         resind, i, j, k;
 783     atom_id     ai, aj;
 784     real        dist2, long_bond_dist2, short_bond_dist2;
 785     const char *ptr;
 786
 787     long_bond_dist2  = sqr(long_bond_dist);
 788     short_bond_dist2 = sqr(short_bond_dist);
 789
 790     if (debug)
 791     {
 792         ptr = "bond";
 793     }
 794     else
 795     {
 796         ptr = "check";
 797     }
 798
 799     fprintf(stderr, "Making bonds...\n");
 800     i = 0;
 801     for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
 802     {
 803         /* add bonds from list of bonded interactions */
 804         for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
 805         {
 806             /* Unfortunately we can not issue errors or warnings
 807              * for missing atoms in bonds, as the hydrogens and terminal atoms
 808              * have not been added yet.
 809              */
 810             ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[0], i, atoms,
 811                              ptr, TRUE);
 812             aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[1], i, atoms,
 813                              ptr, TRUE);
 814             if (ai != NO_ATID && aj != NO_ATID)
 815             {
 816                 dist2 = distance2(x[ai], x[aj]);
 817                 if (dist2 > long_bond_dist2)
 818                 {
 819                     fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
 820                             ai+1, aj+1, sqrt(dist2));
 821                 }
 822                 else if (dist2 < short_bond_dist2)
 823                 {
 824                     fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
 825                             ai+1, aj+1, sqrt(dist2));
 826                 }
 827                 add_param(psb, ai, aj, NULL, hb[resind].rb[ebtsBONDS].b[j].s);
 828             }
 829         }
 830         /* add bonds from list of hacks (each added atom gets a bond) */
 831         while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
 832         {
 833             for (j = 0; j < hb[resind].nhack; j++)
 834             {
 835                 if ( ( hb[resind].hack[j].tp > 0 ||
 836                        hb[resind].hack[j].oname == NULL ) &&
 837                      strcmp(hb[resind].hack[j].a[0], *(atoms->atomname[i])) == 0)
 838                 {
 839                     switch (hb[resind].hack[j].tp)
 840                     {
 841                         case 9:                                   /* COOH terminus */
 842                             add_param(psb, i, i+1, NULL, NULL);   /* C-O  */
 843                             add_param(psb, i, i+2, NULL, NULL);   /* C-OA */
 844                             add_param(psb, i+2, i+3, NULL, NULL); /* OA-H */
 845                             break;
 846                         default:
 847                             for (k = 0; (k < hb[resind].hack[j].nr); k++)
 848                             {
 849                                 add_param(psb, i, i+k+1, NULL, NULL);
 850                             }
 851                     }
 852                 }
 853             }
 854             i++;
 855         }
 856         /* we're now at the start of the next residue */
 857     }
 858 }
 859
 860 static int pcompar(const void *a, const void *b)
 861 {
 862     t_param *pa, *pb;
 863     int      d;
 864     pa = (t_param *)a;
 865     pb = (t_param *)b;
 866
 867     d = pa->a[0] - pb->a[0];
 868     if (d == 0)
 869     {
 870         d = pa->a[1] - pb->a[1];
 871     }
 872     if (d == 0)
 873     {
 874         return strlen(pb->s) - strlen(pa->s);
 875     }
 876     else
 877     {
 878         return d;
 879     }
 880 }
 881
 882 static void clean_bonds(t_params *ps)
 883 {
 884     int     i, j;
 885     atom_id a;
 886
 887     if (ps->nr > 0)
 888     {
 889         /* swap atomnumbers in bond if first larger than second: */
 890         for (i = 0; (i < ps->nr); i++)
 891         {
 892             if (ps->param[i].a[1] < ps->param[i].a[0])
 893             {
 894                 a                 = ps->param[i].a[0];
 895                 ps->param[i].a[0] = ps->param[i].a[1];
 896                 ps->param[i].a[1] = a;
 897             }
 898         }
 899
 900         /* Sort bonds */
 901         qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
 902
 903         /* remove doubles, keep the first one always. */
 904         j = 1;
 905         for (i = 1; (i < ps->nr); i++)
 906         {
 907             if ((ps->param[i].a[0] != ps->param[j-1].a[0]) ||
 908                 (ps->param[i].a[1] != ps->param[j-1].a[1]) )
 909             {
 910                 if (j != i)
 911                 {
 912                     cp_param(&(ps->param[j]), &(ps->param[i]));
 913                 }
 914                 j++;
 915             }
 916         }
 917         fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
 918         ps->nr = j;
 919     }
 920     else
 921     {
 922         fprintf(stderr, "No bonds\n");
 923     }
 924 }
 925
 926 void print_sums(t_atoms *atoms, gmx_bool bSystem)
 927 {
 928     double      m, qtot;
 929     int         i;
 930     const char *where;
 931
 932     if (bSystem)
 933     {
 934         where = " in system";
 935     }
 936     else
 937     {
 938         where = "";
 939     }
 940
 941     m    = 0;
 942     qtot = 0;
 943     for (i = 0; (i < atoms->nr); i++)
 944     {
 945         m    += atoms->atom[i].m;
 946         qtot += atoms->atom[i].q;
 947     }
 948
 949     fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
 950     fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
 951 }
 952
 953 static void check_restp_type(const char *name, int t1, int t2)
 954 {
 955     if (t1 != t2)
 956     {
 957         gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
 958     }
 959 }
 960
 961 static void check_restp_types(t_restp *r0, t_restp *r1)
 962 {
 963     int i;
 964
 965     check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
 966     check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
 967     check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
 968     check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
 969
 970     for (i = 0; i < ebtsNR; i++)
 971     {
 972         check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
 973     }
 974 }
 975
 976 void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
 977 {
 978     char        buf[STRLEN];
 979     int         k;
 980     const char *Hnum = "123456";
 981
 982     /*if (debug)
 983        {
 984         fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
 985                 hack->nname,
 986      * restp->atomname[at_start], resnr, restp->resname);
 987                 }*/
 988     strcpy(buf, hack->nname);
 989     buf[strlen(buf)+1] = '\0';
 990     if (hack->nr > 1)
 991     {
 992         buf[strlen(buf)] = '-';
 993     }
 994     /* make space */
 995     restp->natom += hack->nr;
 996     srenew(restp->atom,     restp->natom);
 997     srenew(restp->atomname, restp->natom);
 998     srenew(restp->cgnr,     restp->natom);
 999     /* shift the rest */
1000     for (k = restp->natom-1; k > at_start+hack->nr; k--)
1001     {
1002         restp->atom[k] =
1003             restp->atom    [k - hack->nr];
1004         restp->atomname[k] =
1005             restp->atomname[k - hack->nr];
1006         restp->cgnr[k] =
1007             restp->cgnr    [k - hack->nr];
1008     }
1009     /* now add them */
1010     for (k = 0; k < hack->nr; k++)
1011     {
1012         /* set counter in atomname */
1013         if (hack->nr > 1)
1014         {
1015             buf[strlen(buf)-1] = Hnum[k];
1016         }
1017         snew( restp->atomname[at_start+1+k], 1);
1018         restp->atom     [at_start+1+k] = *hack->atom;
1019         *restp->atomname[at_start+1+k] = strdup(buf);
1020         if (hack->cgnr != NOTSET)
1021         {
1022             restp->cgnr   [at_start+1+k] = hack->cgnr;
1023         }
1024         else
1025         {
1026             restp->cgnr   [at_start+1+k] = restp->cgnr[at_start];
1027         }
1028     }
1029 }
1030
1031 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1032                         int nrtp, t_restp rtp[],
1033                         int nres, t_resinfo *resinfo,
1034                         int nterpairs,
1035                         t_hackblock **ntdb, t_hackblock **ctdb,
1036                         int *rn, int *rc)
1037 {
1038     int         i, j, k, l;
1039     char       *key;
1040     t_restp    *res;
1041     int         tern, terc;
1042     gmx_bool    bRM;
1043
1044     snew(*hb, nres);
1045     snew(*restp, nres);
1046     /* first the termini */
1047     for (i = 0; i < nterpairs; i++)
1048     {
1049         if (rn[i] >= 0 && ntdb[i] != NULL)
1050         {
1051             copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1052         }
1053         if (rc[i] >= 0 && ctdb[i] != NULL)
1054         {
1055             merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1056         }
1057     }
1058
1059     /* then the whole rtp */
1060     for (i = 0; i < nres; i++)
1061     {
1062         /* Here we allow a mismatch of one character when looking for the rtp entry.
1063          * For such a mismatch there should be only one mismatching name.
1064          * This is mainly useful for small molecules such as ions.
1065          * Note that this will usually not work for protein, DNA and RNA,
1066          * since there the residue names should be listed in residuetypes.dat
1067          * and an error will have been generated earlier in the process.
1068          */
1069         key = *resinfo[i].rtp;
1070         snew(resinfo[i].rtp, 1);
1071         *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1072         res             = get_restp(*resinfo[i].rtp, nrtp, rtp);
1073         copy_t_restp(res, &(*restp)[i]);
1074
1075         /* Check that we do not have different bonded types in one molecule */
1076         check_restp_types(&(*restp)[0], &(*restp)[i]);
1077
1078         tern = -1;
1079         for (j = 0; j < nterpairs && tern == -1; j++)
1080         {
1081             if (i == rn[j])
1082             {
1083                 tern = j;
1084             }
1085         }
1086         terc = -1;
1087         for (j = 0; j < nterpairs && terc == -1; j++)
1088         {
1089             if (i == rc[j])
1090             {
1091                 terc = j;
1092             }
1093         }
1094         bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1095
1096         if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1097                     (terc >= 0 && ctdb[terc] == NULL)))
1098         {
1099             gmx_fatal(FARGS, "There is a dangling bond at at least one of the terminal ends and the force field does not provide terminal entries or files. Fix your terminal residues so that they match the residue database (.rtp) entries, or provide terminal database entries (.tdb).");
1100         }
1101         if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1102                     (terc >= 0 && ctdb[terc]->nhack == 0)))
1103         {
1104             gmx_fatal(FARGS, "There is a dangling bond at at least one of the terminal ends. Fix your coordinate file, add a new terminal database entry (.tdb), or select the proper existing terminal entry.");
1105         }
1106     }
1107
1108     /* now perform t_hack's on t_restp's,
1109        i.e. add's and deletes from termini database will be
1110        added to/removed from residue topology
1111        we'll do this on one big dirty loop, so it won't make easy reading! */
1112     for (i = 0; i < nres; i++)
1113     {
1114         for (j = 0; j < (*hb)[i].nhack; j++)
1115         {
1116             if ( (*hb)[i].hack[j].nr)
1117             {
1118                 /* find atom in restp */
1119                 for (l = 0; l < (*restp)[i].natom; l++)
1120                 {
1121                     if ( ( (*hb)[i].hack[j].oname == NULL &&
1122                            strcmp((*hb)[i].hack[j].a[0], *(*restp)[i].atomname[l]) == 0 ) ||
1123                          ( (*hb)[i].hack[j].oname != NULL &&
1124                            strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1125                     {
1126                         break;
1127                     }
1128                 }
1129                 if (l == (*restp)[i].natom)
1130                 {
1131                     /* If we are doing an atom rename only, we don't need
1132                      * to generate a fatal error if the old name is not found
1133                      * in the rtp.
1134                      */
1135                     /* Deleting can happen also only on the input atoms,
1136                      * not necessarily always on the rtp entry.
1137                      */
1138                     if (!((*hb)[i].hack[j].oname != NULL &&
1139                           (*hb)[i].hack[j].nname != NULL) &&
1140                         !((*hb)[i].hack[j].oname != NULL &&
1141                           (*hb)[i].hack[j].nname == NULL))
1142                     {
1143                         gmx_fatal(FARGS,
1144                                   "atom %s not found in buiding block %d%s "
1145                                   "while combining tdb and rtp",
1146                                   (*hb)[i].hack[j].oname != NULL ?
1147                                   (*hb)[i].hack[j].oname : (*hb)[i].hack[j].a[0],
1148                                   i+1, *resinfo[i].rtp);
1149                     }
1150                 }
1151                 else
1152                 {
1153                     if ( (*hb)[i].hack[j].oname == NULL)
1154                     {
1155                         /* we're adding: */
1156                         add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1157                     }
1158                     else
1159                     {
1160                         /* oname != NULL */
1161                         if ( (*hb)[i].hack[j].nname == NULL)
1162                         {
1163                             /* we're deleting */
1164                             if (debug)
1165                             {
1166                                 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1167                                         *(*restp)[i].atomname[l],
1168                                         i+1, (*restp)[i].resname);
1169                             }
1170                             /* shift the rest */
1171                             (*restp)[i].natom--;
1172                             for (k = l; k < (*restp)[i].natom; k++)
1173                             {
1174                                 (*restp)[i].atom    [k] = (*restp)[i].atom    [k+1];
1175                                 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1176                                 (*restp)[i].cgnr    [k] = (*restp)[i].cgnr    [k+1];
1177                             }
1178                             /* give back space */
1179                             srenew((*restp)[i].atom,     (*restp)[i].natom);
1180                             srenew((*restp)[i].atomname, (*restp)[i].natom);
1181                             srenew((*restp)[i].cgnr,     (*restp)[i].natom);
1182                         }
1183                         else /* nname != NULL */
1184                         {    /* we're replacing */
1185                             if (debug)
1186                             {
1187                                 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1188                                         *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1189                                         i+1, (*restp)[i].resname);
1190                             }
1191                             snew( (*restp)[i].atomname[l], 1);
1192                             (*restp)[i].atom[l]      =       *(*hb)[i].hack[j].atom;
1193                             *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1194                             if ( (*hb)[i].hack[j].cgnr != NOTSET)
1195                             {
1196                                 (*restp)[i].cgnr   [l] = (*hb)[i].hack[j].cgnr;
1197                             }
1198                         }
1199                     }
1200                 }
1201             }
1202         }
1203     }
1204 }
1205
1206 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1207 {
1208
1209     if (hack->nr == 1)
1210     {
1211         *nr = 0;
1212
1213         return (gmx_strcasecmp(anm, hack->nname) == 0);
1214     }
1215     else
1216     {
1217         if (isdigit(anm[strlen(anm)-1]))
1218         {
1219             *nr = anm[strlen(anm)-1] - '0';
1220         }
1221         else
1222         {
1223             *nr = 0;
1224         }
1225         if (*nr <= 0 || *nr > hack->nr)
1226         {
1227             return FALSE;
1228         }
1229         else
1230         {
1231             return (strlen(anm) == strlen(hack->nname) + 1 &&
1232                     gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1233         }
1234     }
1235 }
1236
1237 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1238                                               t_restp *rptr, t_hackblock *hbr,
1239                                               gmx_bool bVerbose)
1240 {
1241     int      resnr;
1242     int      j, k;
1243     char    *oldnm, *newnm;
1244     int      anmnr;
1245     char    *start_at, buf[STRLEN];
1246     int      start_nr;
1247     gmx_bool bReplaceReplace, bFoundInAdd;
1248     gmx_bool bDeleted;
1249
1250     oldnm = *pdba->atomname[atind];
1251     resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1252
1253     bDeleted = FALSE;
1254     for (j = 0; j < hbr->nhack; j++)
1255     {
1256         if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1257             gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1258         {
1259             /* This is a replace entry. */
1260             /* Check if we are not replacing a replaced atom. */
1261             bReplaceReplace = FALSE;
1262             for (k = 0; k < hbr->nhack; k++)
1263             {
1264                 if (k != j &&
1265                     hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1266                     gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1267                 {
1268                     /* The replace in hack[j] replaces an atom that
1269                      * was already replaced in hack[k], we do not want
1270                      * second or higher level replaces at this stage.
1271                      */
1272                     bReplaceReplace = TRUE;
1273                 }
1274             }
1275             if (bReplaceReplace)
1276             {
1277                 /* Skip this replace. */
1278                 continue;
1279             }
1280
1281             /* This atom still has the old name, rename it */
1282             newnm = hbr->hack[j].nname;
1283             for (k = 0; k < rptr->natom; k++)
1284             {
1285                 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1286                 {
1287                     break;
1288                 }
1289             }
1290             if (k == rptr->natom)
1291             {
1292                 /* The new name is not present in the rtp.
1293                  * We need to apply the replace also to the rtp entry.
1294                  */
1295
1296                 /* We need to find the add hack that can add this atom
1297                  * to find out after which atom it should be added.
1298                  */
1299                 bFoundInAdd = FALSE;
1300                 for (k = 0; k < hbr->nhack; k++)
1301                 {
1302                     if (hbr->hack[k].oname == NULL &&
1303                         hbr->hack[k].nname != NULL &&
1304                         atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1305                     {
1306                         if (anmnr <= 1)
1307                         {
1308                             start_at = hbr->hack[k].a[0];
1309                         }
1310                         else
1311                         {
1312                             sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1313                             start_at = buf;
1314                         }
1315                         for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1316                         {
1317                             if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1318                             {
1319                                 break;
1320                             }
1321                         }
1322                         if (start_nr == rptr->natom)
1323                         {
1324                             gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1325                                       start_at, rptr->resname, newnm);
1326                         }
1327                         /* We can add the atom after atom start_nr */
1328                         add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1329
1330                         bFoundInAdd = TRUE;
1331                     }
1332                 }
1333
1334                 if (!bFoundInAdd)
1335                 {
1336                     gmx_fatal(FARGS, "Could not find an 'add' entry for atom named '%s' corresponding to the 'replace' entry from atom name '%s' to '%s' for tdb or hdb database of residue type '%s'",
1337                               newnm,
1338                               hbr->hack[j].oname, hbr->hack[j].nname,
1339                               rptr->resname);
1340                 }
1341             }
1342
1343             if (bVerbose)
1344             {
1345                 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1346                        oldnm, rptr->resname, resnr, newnm);
1347             }
1348             /* Rename the atom in pdba */
1349             snew(pdba->atomname[atind], 1);
1350             *pdba->atomname[atind] = strdup(newnm);
1351         }
1352         else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1353                  gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1354         {
1355             /* This is a delete entry, check if this atom is present
1356              * in the rtp entry of this residue.
1357              */
1358             for (k = 0; k < rptr->natom; k++)
1359             {
1360                 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1361                 {
1362                     break;
1363                 }
1364             }
1365             if (k == rptr->natom)
1366             {
1367                 /* This atom is not present in the rtp entry,
1368                  * delete is now from the input pdba.
1369                  */
1370                 if (bVerbose)
1371                 {
1372                     printf("Deleting atom '%s' in residue '%s' %d\n",
1373                            oldnm, rptr->resname, resnr);
1374                 }
1375                 /* We should free the atom name,
1376                  * but it might be used multiple times in the symtab.
1377                  * sfree(pdba->atomname[atind]);
1378                  */
1379                 for (k = atind+1; k < pdba->nr; k++)
1380                 {
1381                     pdba->atom[k-1]     = pdba->atom[k];
1382                     pdba->atomname[k-1] = pdba->atomname[k];
1383                     copy_rvec(x[k], x[k-1]);
1384                 }
1385                 pdba->nr--;
1386                 bDeleted = TRUE;
1387             }
1388         }
1389     }
1390
1391     return bDeleted;
1392 }
1393
1394 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1395                               t_atoms *pdba, rvec *x,
1396                               gmx_bool bVerbose)
1397 {
1398     int          i, j;
1399     char        *oldnm;
1400     t_restp     *rptr;
1401
1402     for (i = 0; i < pdba->nr; i++)
1403     {
1404         oldnm = *pdba->atomname[i];
1405         rptr  = &restp[pdba->atom[i].resind];
1406         for (j = 0; (j < rptr->natom); j++)
1407         {
1408             if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1409             {
1410                 break;
1411             }
1412         }
1413         if (j == rptr->natom)
1414         {
1415             /* Not found yet, check if we have to rename this atom */
1416             if (match_atomnames_with_rtp_atom(pdba, x, i,
1417                                               rptr, &(hb[pdba->atom[i].resind]),
1418                                               bVerbose))
1419             {
1420                 /* We deleted this atom, decrease the atom counter by 1. */
1421                 i--;
1422             }
1423         }
1424     }
1425 }
1426
1427 #define NUM_CMAP_ATOMS 5
1428 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t *rt)
1429 {
1430     int         residx, i, j, k;
1431     const char *ptr;
1432     t_resinfo  *resinfo = atoms->resinfo;
1433     int         nres    = atoms->nres;
1434     gmx_bool    bAddCMAP;
1435     atom_id     cmap_atomid[NUM_CMAP_ATOMS];
1436     int         cmap_chainnum = -1, this_residue_index;
1437
1438     if (debug)
1439     {
1440         ptr = "cmap";
1441     }
1442     else
1443     {
1444         ptr = "check";
1445     }
1446
1447     fprintf(stderr, "Making cmap torsions...");
1448     i = 0;
1449     /* End loop at nres-1, since the very last residue does not have a +N atom, and
1450      * therefore we get a valgrind invalid 4 byte read error with atom am */
1451     for (residx = 0; residx < nres-1; residx++)
1452     {
1453         /* Add CMAP terms from the list of CMAP interactions */
1454         for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1455         {
1456             bAddCMAP = TRUE;
1457             /* Loop over atoms in a candidate CMAP interaction and
1458              * check that they exist, are from the same chain and are
1459              * from residues labelled as protein. */
1460             for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1461             {
1462                 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
1463                                              i, atoms, ptr, TRUE);
1464                 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1465                 if (!bAddCMAP)
1466                 {
1467                     /* This break is necessary, because cmap_atomid[k]
1468                      * == NO_ATID cannot be safely used as an index
1469                      * into the atom array. */
1470                     break;
1471                 }
1472                 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1473                 if (0 == k)
1474                 {
1475                     cmap_chainnum = resinfo[this_residue_index].chainnum;
1476                 }
1477                 else
1478                 {
1479                     /* Does the residue for this atom have the same
1480                      * chain number as the residues for previous
1481                      * atoms? */
1482                     bAddCMAP = bAddCMAP &&
1483                         cmap_chainnum == resinfo[this_residue_index].chainnum;
1484                 }
1485                 bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt, *(resinfo[this_residue_index].name));
1486             }
1487
1488             if (bAddCMAP)
1489             {
1490                 add_cmap_param(psb, cmap_atomid[0], cmap_atomid[1], cmap_atomid[2], cmap_atomid[3], cmap_atomid[4], restp[residx].rb[ebtsCMAP].b[j].s);
1491             }
1492         }
1493
1494         if (residx < nres-1)
1495         {
1496             while (atoms->atom[i].resind < residx+1)
1497             {
1498                 i++;
1499             }
1500         }
1501     }
1502
1503     /* Start the next residue */
1504 }
1505
1506 static void
1507 scrub_charge_groups(int *cgnr, int natoms)
1508 {
1509     int i;
1510
1511     for (i = 0; i < natoms; i++)
1512     {
1513         cgnr[i] = i+1;
1514     }
1515 }
1516
1517
1518 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1519              t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1520              int nrtp, t_restp rtp[],
1521              t_restp *restp, t_hackblock *hb,
1522              gmx_bool bAllowMissing,
1523              gmx_bool bVsites, gmx_bool bVsiteAromatics,
1524              const char *ffdir,
1525              real mHmult,
1526              int nssbonds, t_ssbond *ssbonds,
1527              real long_bond_dist, real short_bond_dist,
1528              gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1529              gmx_bool bRenumRes, gmx_bool bRTPresname)
1530 {
1531     /*
1532        t_hackblock *hb;
1533        t_restp  *restp;
1534      */
1535     t_params          plist[F_NRE];
1536     t_excls          *excls;
1537     t_nextnb          nnb;
1538     int              *cgnr;
1539     int              *vsite_type;
1540     int               i, nmissat;
1541     int               bts[ebtsNR];
1542     gmx_residuetype_t*rt;
1543
1544     init_plist(plist);
1545     gmx_residuetype_init(&rt);
1546
1547     if (debug)
1548     {
1549         print_resall(debug, atoms->nres, restp, atype);
1550         dump_hb(debug, atoms->nres, hb);
1551     }
1552
1553     /* Make bonds */
1554     at2bonds(&(plist[F_BONDS]), hb,
1555              atoms, *x,
1556              long_bond_dist, short_bond_dist);
1557
1558     /* specbonds: disulphide bonds & heme-his */
1559     do_ssbonds(&(plist[F_BONDS]),
1560                atoms, nssbonds, ssbonds,
1561                bAllowMissing);
1562
1563     nmissat = name2type(atoms, &cgnr, restp, rt);
1564     if (nmissat)
1565     {
1566         if (bAllowMissing)
1567         {
1568             fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1569                     nmissat, molname);
1570         }
1571         else
1572         {
1573             gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1574                       nmissat, molname);
1575         }
1576     }
1577
1578     /* Cleanup bonds (sort and rm doubles) */
1579     clean_bonds(&(plist[F_BONDS]));
1580
1581     snew(vsite_type, atoms->nr);
1582     for (i = 0; i < atoms->nr; i++)
1583     {
1584         vsite_type[i] = NOTSET;
1585     }
1586     if (bVsites)
1587     {
1588         /* determine which atoms will be vsites and add dummy masses
1589            also renumber atom numbers in plist[0..F_NRE]! */
1590         do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1591                   &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1592     }
1593
1594     /* Make Angles and Dihedrals */
1595     fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1596     snew(excls, atoms->nr);
1597     init_nnb(&nnb, atoms->nr, 4);
1598     gen_nnb(&nnb, plist);
1599     print_nnb(&nnb, "NNB");
1600     gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1601     done_nnb(&nnb);
1602
1603     /* Make CMAP */
1604     if (TRUE == bCmap)
1605     {
1606         gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1607         if (plist[F_CMAP].nr > 0)
1608         {
1609             fprintf(stderr, "There are %4d cmap torsion pairs\n",
1610                     plist[F_CMAP].nr);
1611         }
1612     }
1613
1614     /* set mass of all remaining hydrogen atoms */
1615     if (mHmult != 1.0)
1616     {
1617         do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1618     }
1619     sfree(vsite_type);
1620
1621     /* Cleanup bonds (sort and rm doubles) */
1622     /* clean_bonds(&(plist[F_BONDS]));*/
1623
1624     fprintf(stderr,
1625             "There are %4d dihedrals, %4d impropers, %4d angles\n"
1626             "          %4d pairs,     %4d bonds and  %4d virtual sites\n",
1627             plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1628             plist[F_LJ14].nr, plist[F_BONDS].nr,
1629             plist[F_VSITE2].nr +
1630             plist[F_VSITE3].nr +
1631             plist[F_VSITE3FD].nr +
1632             plist[F_VSITE3FAD].nr +
1633             plist[F_VSITE3OUT].nr +
1634             plist[F_VSITE4FD].nr +
1635             plist[F_VSITE4FDN].nr );
1636
1637     print_sums(atoms, FALSE);
1638
1639     if (FALSE == bChargeGroups)
1640     {
1641         scrub_charge_groups(cgnr, atoms->nr);
1642     }
1643
1644     if (bRenumRes)
1645     {
1646         for (i = 0; i < atoms->nres; i++)
1647         {
1648             atoms->resinfo[i].nr = i + 1;
1649             atoms->resinfo[i].ic = ' ';
1650         }
1651     }
1652
1653     if (top_file)
1654     {
1655         fprintf(stderr, "Writing topology\n");
1656         /* We can copy the bonded types from the first restp,
1657          * since the types have to be identical for all residues in one molecule.
1658          */
1659         for (i = 0; i < ebtsNR; i++)
1660         {
1661             bts[i] = restp[0].rb[i].type;
1662         }
1663         write_top(top_file, posre_fn, molname,
1664                   atoms, bRTPresname,
1665                   bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1666     }
1667
1668     /* cleaning up */
1669     free_t_hackblock(atoms->nres, &hb);
1670     free_t_restp(atoms->nres, &restp);
1671     gmx_residuetype_destroy(rt);
1672
1673     /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1674     sfree(cgnr);
1675     for (i = 0; i < F_NRE; i++)
1676     {
1677         sfree(plist[i].param);
1678     }
1679     for (i = 0; i < atoms->nr; i++)
1680     {
1681         sfree(excls[i].e);
1682     }
1683     sfree(excls);
1684 }