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