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