src/gromacs/gmxpreprocess/pdb2top.cpp

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