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