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