src/gromacs/gmxpreprocess/pdb2top.cpp

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