src/gromacs/gmxpreprocess/pdb2top.cpp

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