src/gromacs/gmxpreprocess/pdb2top.cpp

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