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