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