src/gromacs/gmxpreprocess/pdb2top.cpp

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