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