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