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