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