src/gromacs/gmxpreprocess/pdb2top.c

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