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