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