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                PreprocessingAtomTypes *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);
 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, {}, 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);   /* C-O  */
 799                             add_param(psb, i, i+2, {}, nullptr);   /* C-OA */
 800                             add_param(psb, i+2, i+3, {}, 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);
 806                             }
 807                     }
 808                 }
 809             }
 810             i++;
 811         }
 812         /* we're now at the start of the next residue */
 813     }
 814 }
 815
 816 static bool pcompar(const InteractionType &a, const InteractionType &b)
 817 {
 818     int                d;
 819
 820     if ((d = a.ai() - b.ai()) != 0)
 821     {
 822         return d < 0;
 823     }
 824     else if ((d = a.aj() - b.aj()) != 0)
 825     {
 826         return d < 0;
 827     }
 828     else
 829     {
 830         return a.interactionTypeName().length() > b.interactionTypeName().length();
 831     }
 832 }
 833
 834 static void clean_bonds(InteractionTypeParameters *ps)
 835 {
 836     if (ps->size() > 0)
 837     {
 838         /* Sort bonds */
 839         for (auto &bond : ps->interactionTypes)
 840         {
 841             bond.sortAtomIds();
 842         }
 843         std::sort(ps->interactionTypes.begin(), ps->interactionTypes.end(), pcompar);
 844
 845         /* remove doubles, keep the first one always. */
 846         int oldNumber = ps->size();
 847         for (auto parm = ps->interactionTypes.begin() + 1; parm != ps->interactionTypes.end(); )
 848         {
 849             auto prev = parm - 1;
 850             if (parm->ai() == prev->ai() &&
 851                 parm->aj() == prev->aj())
 852             {
 853                 parm = ps->interactionTypes.erase(parm);
 854             }
 855             else
 856             {
 857                 ++parm;
 858             }
 859         }
 860         fprintf(stderr, "Number of bonds was %d, now %zu\n", oldNumber, ps->size());
 861     }
 862     else
 863     {
 864         fprintf(stderr, "No bonds\n");
 865     }
 866 }
 867
 868 void print_sums(const t_atoms *atoms, bool bSystem)
 869 {
 870     double      m, qtot;
 871     int         i;
 872     const char *where;
 873
 874     if (bSystem)
 875     {
 876         where = " in system";
 877     }
 878     else
 879     {
 880         where = "";
 881     }
 882
 883     m    = 0;
 884     qtot = 0;
 885     for (i = 0; (i < atoms->nr); i++)
 886     {
 887         m    += atoms->atom[i].m;
 888         qtot += atoms->atom[i].q;
 889     }
 890
 891     fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
 892     fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
 893 }
 894
 895 static void check_restp_type(const char *name, int t1, int t2)
 896 {
 897     if (t1 != t2)
 898     {
 899         gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
 900     }
 901 }
 902
 903 static void check_restp_types(const PreprocessResidue &r0, const PreprocessResidue &r1)
 904 {
 905     check_restp_type("all dihedrals", static_cast<int>(r0.bKeepAllGeneratedDihedrals), static_cast<int>(r1.bKeepAllGeneratedDihedrals));
 906     check_restp_type("nrexcl", r0.nrexcl, r1.nrexcl);
 907     check_restp_type("HH14", static_cast<int>(r0.bGenerateHH14Interactions), static_cast<int>(r1.bGenerateHH14Interactions));
 908     check_restp_type("remove dihedrals", static_cast<int>(r0.bRemoveDihedralIfWithImproper), static_cast<int>(r1.bRemoveDihedralIfWithImproper));
 909
 910     for (int i = 0; i < ebtsNR; i++)
 911     {
 912         check_restp_type(btsNames[i], r0.rb[i].type, r1.rb[i].type);
 913     }
 914 }
 915
 916 static void add_atom_to_restp(PreprocessResidue         *usedPpResidues,
 917                               t_symtab                  *symtab,
 918                               int                        at_start,
 919                               const MoleculePatch       *patch)
 920 {
 921     /* now add them */
 922     for (int k = 0; k < patch->nr; k++)
 923     {
 924         /* set counter in atomname */
 925         std::string buf = patch->nname;
 926         if (patch->nr > 1)
 927         {
 928             buf.append(gmx::formatString("%d", k+1));
 929         }
 930         usedPpResidues->atomname.insert(usedPpResidues->atomname.begin() + at_start + 1 + k,
 931                                         put_symtab(symtab, buf.c_str()));
 932         usedPpResidues->atom.insert(usedPpResidues->atom.begin() + at_start + 1 + k, patch->atom.back());
 933         if (patch->cgnr != NOTSET)
 934         {
 935             usedPpResidues->cgnr.insert(usedPpResidues->cgnr.begin() + at_start + 1 + k, patch->cgnr);
 936         }
 937         else
 938         {
 939             usedPpResidues->cgnr.insert(usedPpResidues->cgnr.begin() + at_start + 1 + k, usedPpResidues->cgnr[at_start]);
 940         }
 941     }
 942 }
 943
 944 void get_hackblocks_rtp(std::vector<MoleculePatchDatabase> *globalPatches,
 945                         std::vector<PreprocessResidue> *usedPpResidues,
 946                         gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
 947                         int nres, t_resinfo *resinfo,
 948                         int nterpairs,
 949                         t_symtab *symtab,
 950                         gmx::ArrayRef<MoleculePatchDatabase *> ntdb,
 951                         gmx::ArrayRef<MoleculePatchDatabase *> ctdb,
 952                         gmx::ArrayRef<const int> rn,
 953                         gmx::ArrayRef<const int> rc,
 954                         bool bAllowMissing)
 955 {
 956     char       *key;
 957     bool        bRM;
 958
 959     globalPatches->resize(nres);
 960     usedPpResidues->clear();
 961     /* first the termini */
 962     for (int i = 0; i < nterpairs; i++)
 963     {
 964         if (rn[i] >= 0 && ntdb[i] != nullptr)
 965         {
 966             copyModificationBlocks(*ntdb[i], &globalPatches->at(rn[i]));
 967         }
 968         if (rc[i] >= 0 && ctdb[i] != nullptr)
 969         {
 970             mergeAtomAndBondModifications(*ctdb[i], &globalPatches->at(rc[i]));
 971         }
 972     }
 973
 974     /* then the whole rtp */
 975     for (int i = 0; i < nres; i++)
 976     {
 977         /* Here we allow a mismatch of one character when looking for the rtp entry.
 978          * For such a mismatch there should be only one mismatching name.
 979          * This is mainly useful for small molecules such as ions.
 980          * Note that this will usually not work for protein, DNA and RNA,
 981          * since there the residue names should be listed in residuetypes.dat
 982          * and an error will have been generated earlier in the process.
 983          */
 984         key = *resinfo[i].rtp;
 985
 986         resinfo[i].rtp = put_symtab(symtab, searchResidueDatabase(key, rtpFFDB).c_str());
 987         auto               res             = getDatabaseEntry(*resinfo[i].rtp, rtpFFDB);
 988         usedPpResidues->push_back(PreprocessResidue());
 989         PreprocessResidue *newentry = &usedPpResidues->back();
 990         copyPreprocessResidues(*res, newentry, symtab);
 991
 992         /* Check that we do not have different bonded types in one molecule */
 993         check_restp_types(usedPpResidues->front(), *newentry);
 994
 995         int tern = -1;
 996         for (int j = 0; j < nterpairs && tern == -1; j++)
 997         {
 998             if (i == rn[j])
 999             {
1000                 tern = j;
1001             }
1002         }
1003         int terc = -1;
1004         for (int j = 0; j < nterpairs && terc == -1; j++)
1005         {
1006             if (i == rc[j])
1007             {
1008                 terc = j;
1009             }
1010         }
1011         bRM = mergeBondedInteractionList(res->rb, globalPatches->at(i).rb, tern >= 0, terc >= 0);
1012
1013         if (bRM && ((tern >= 0 && ntdb[tern] == nullptr) ||
1014                     (terc >= 0 && ctdb[terc] == nullptr)))
1015         {
1016             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).";
1017             if (bAllowMissing)
1018             {
1019                 fprintf(stderr, "%s\n", errString);
1020             }
1021             else
1022             {
1023                 gmx_fatal(FARGS, "%s", errString);
1024             }
1025         }
1026         else if (bRM && ((tern >= 0 && ntdb[tern]->nhack() == 0) ||
1027                          (terc >= 0 && ctdb[terc]->nhack() == 0)))
1028         {
1029             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.";
1030             if (bAllowMissing)
1031             {
1032                 fprintf(stderr, "%s\n", errString);
1033             }
1034             else
1035             {
1036                 gmx_fatal(FARGS, "%s", errString);
1037             }
1038         }
1039     }
1040
1041     /* Apply patchs to t_restp entries
1042        i.e. add's and deletes from termini database will be
1043        added to/removed from residue topology
1044        we'll do this on one big dirty loop, so it won't make easy reading! */
1045     for (auto modifiedResidue = globalPatches->begin();
1046          modifiedResidue != globalPatches->end();
1047          modifiedResidue++)
1048     {
1049         const int pos = std::distance(globalPatches->begin(), modifiedResidue);
1050         PreprocessResidue             *posres = &usedPpResidues->at(pos);
1051         for (auto patch = modifiedResidue->hack.begin();
1052              patch != modifiedResidue->hack.end();
1053              patch++)
1054         {
1055             if (patch->nr != 0)
1056             {
1057                 /* find atom in restp */
1058                 auto found =
1059                     std::find_if(posres->atomname.begin(),
1060                                  posres->atomname.end(),
1061                                  [&patch](char **name)
1062                                  { return (patch->oname.empty() &&
1063                                            patch->a[0] == *name) ||
1064                                    (patch->oname == *name); });
1065
1066                 if (found == posres->atomname.end())
1067                 {
1068                     /* If we are doing an atom rename only, we don't need
1069                      * to generate a fatal error if the old name is not found
1070                      * in the rtp.
1071                      */
1072                     /* Deleting can happen also only on the input atoms,
1073                      * not necessarily always on the rtp entry.
1074                      */
1075                     if (patch->type() == MoleculePatchType::Add)
1076                     {
1077                         gmx_fatal(FARGS,
1078                                   "atom %s not found in buiding block %d%s "
1079                                   "while combining tdb and rtp",
1080                                   patch->oname.empty() ?
1081                                   patch->a[0].c_str() : patch->oname.c_str(),
1082                                   pos+1, *resinfo[pos].rtp);
1083                     }
1084                 }
1085                 else
1086                 {
1087                     int l = std::distance(posres->atomname.begin(), found);
1088                     switch (patch->type())
1089                     {
1090                         case MoleculePatchType::Add:
1091                         {
1092                             /* we're adding: */
1093                             add_atom_to_restp(posres, symtab, l, &(*patch));
1094                             break;
1095                         }
1096                         case MoleculePatchType::Delete:
1097                         {   /* we're deleting */
1098                             posres->atom.erase(posres->atom.begin() + l);
1099                             posres->atomname.erase(posres->atomname.begin() + l);
1100                             posres->cgnr.erase(posres->cgnr.begin() + l);
1101                             break;
1102                         }
1103                         case MoleculePatchType::Replace:
1104                         {
1105                             /* we're replacing */
1106                             posres->atom[l]      = patch->atom.back();
1107                             posres->atomname[l]  = put_symtab(symtab, patch->nname.c_str());
1108                             if (patch->cgnr != NOTSET)
1109                             {
1110                                 posres->cgnr[l] = patch->cgnr;
1111                             }
1112                             break;
1113                         }
1114                     }
1115                 }
1116             }
1117         }
1118     }
1119 }
1120
1121 static bool atomname_cmp_nr(const char *anm, const MoleculePatch *patch, int *nr)
1122 {
1123
1124     if (patch->nr == 1)
1125     {
1126         *nr = 0;
1127
1128         return (gmx::equalCaseInsensitive(anm, patch->nname));
1129     }
1130     else
1131     {
1132         if (isdigit(anm[strlen(anm)-1]))
1133         {
1134             *nr = anm[strlen(anm)-1] - '0';
1135         }
1136         else
1137         {
1138             *nr = 0;
1139         }
1140         if (*nr <= 0 || *nr > patch->nr)
1141         {
1142             return FALSE;
1143         }
1144         else
1145         {
1146             return (strlen(anm) == patch->nname.length() + 1 &&
1147                     gmx_strncasecmp(anm, patch->nname.c_str(), patch->nname.length()) == 0);
1148         }
1149     }
1150 }
1151
1152 static bool match_atomnames_with_rtp_atom(t_atoms                     *pdba,
1153                                           gmx::ArrayRef<gmx::RVec>     x,
1154                                           t_symtab                    *symtab,
1155                                           int                          atind,
1156                                           PreprocessResidue           *localPpResidue,
1157                                           const MoleculePatchDatabase &singlePatch,
1158                                           bool                         bVerbose)
1159 {
1160     int      resnr;
1161     char    *oldnm;
1162     int      anmnr;
1163     bool     bDeleted;
1164
1165     oldnm = *pdba->atomname[atind];
1166     resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1167
1168     bDeleted = FALSE;
1169     for (auto patch = singlePatch.hack.begin();
1170          patch != singlePatch.hack.end();
1171          patch++)
1172     {
1173         if (patch->type() == MoleculePatchType::Replace &&
1174             gmx::equalCaseInsensitive(oldnm, patch->oname))
1175         {
1176             /* This is a replace entry. */
1177             /* Check if we are not replacing a replaced atom. */
1178             bool bReplaceReplace = false;
1179             for (auto selfPatch = singlePatch.hack.begin();
1180                  selfPatch != singlePatch.hack.end();
1181                  selfPatch++)
1182             {
1183                 if (patch != selfPatch &&
1184                     selfPatch->type() == MoleculePatchType::Replace &&
1185                     gmx::equalCaseInsensitive(selfPatch->nname, patch->oname))
1186                 {
1187                     /* The replace in patch replaces an atom that
1188                      * was already replaced in selfPatch, we do not want
1189                      * second or higher level replaces at this stage.
1190                      */
1191                     bReplaceReplace = true;
1192                 }
1193             }
1194             if (bReplaceReplace)
1195             {
1196                 /* Skip this replace. */
1197                 continue;
1198             }
1199
1200             /* This atom still has the old name, rename it */
1201             std::string newnm = patch->nname;
1202             auto        found = std::find_if(localPpResidue->atomname.begin(),
1203                                              localPpResidue->atomname.end(),
1204                                              [&newnm](char** name)
1205                                              { return gmx::equalCaseInsensitive(newnm, *name); });
1206             if (found == localPpResidue->atomname.end())
1207             {
1208                 /* The new name is not present in the rtp.
1209                  * We need to apply the replace also to the rtp entry.
1210                  */
1211
1212                 /* We need to find the add hack that can add this atom
1213                  * to find out after which atom it should be added.
1214                  */
1215                 bool bFoundInAdd = false;
1216                 for (auto rtpModification = singlePatch.hack.begin();
1217                      rtpModification != singlePatch.hack.end();
1218                      rtpModification++)
1219                 {
1220                     int         k = std::distance(localPpResidue->atomname.begin(), found);
1221                     std::string start_at;
1222                     if (rtpModification->type() == MoleculePatchType::Add &&
1223                         atomname_cmp_nr(newnm.c_str(), &(*rtpModification), &anmnr))
1224                     {
1225                         if (anmnr <= 1)
1226                         {
1227                             start_at = singlePatch.hack[k].a[0];
1228                         }
1229                         else
1230                         {
1231                             start_at = gmx::formatString("%s%d", singlePatch.hack[k].nname.c_str(), anmnr-1);
1232                         }
1233                         auto found2 = std::find_if(localPpResidue->atomname.begin(),
1234                                                    localPpResidue->atomname.end(),
1235                                                    [&start_at](char **name)
1236                                                    { return gmx::equalCaseInsensitive(start_at, *name); });
1237                         if (found2 == localPpResidue->atomname.end())
1238                         {
1239                             gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1240                                       start_at.c_str(), localPpResidue->resname.c_str(), newnm.c_str());
1241                         }
1242                         /* We can add the atom after atom start_nr */
1243                         add_atom_to_restp(localPpResidue, symtab,
1244                                           std::distance(localPpResidue->atomname.begin(), found2),
1245                                           &(*patch));
1246
1247                         bFoundInAdd = true;
1248                     }
1249                 }
1250
1251                 if (!bFoundInAdd)
1252                 {
1253                     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'",
1254                               newnm.c_str(),
1255                               patch->oname.c_str(), patch->nname.c_str(),
1256                               localPpResidue->resname.c_str());
1257                 }
1258             }
1259
1260             if (bVerbose)
1261             {
1262                 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1263                        oldnm, localPpResidue->resname.c_str(), resnr, newnm.c_str());
1264             }
1265             /* Rename the atom in pdba */
1266             pdba->atomname[atind] = put_symtab(symtab, newnm.c_str());
1267         }
1268         else if (patch->type() == MoleculePatchType::Delete &&
1269                  gmx::equalCaseInsensitive(oldnm, patch->oname))
1270         {
1271             /* This is a delete entry, check if this atom is present
1272              * in the rtp entry of this residue.
1273              */
1274             auto found3 = std::find_if(localPpResidue->atomname.begin(), localPpResidue->atomname.end(),
1275                                        [&oldnm](char **name)
1276                                        { return gmx::equalCaseInsensitive(oldnm, *name); });
1277             if (found3 == localPpResidue->atomname.end())
1278             {
1279                 /* This atom is not present in the rtp entry,
1280                  * delete is now from the input pdba.
1281                  */
1282                 if (bVerbose)
1283                 {
1284                     printf("Deleting atom '%s' in residue '%s' %d\n",
1285                            oldnm, localPpResidue->resname.c_str(), resnr);
1286                 }
1287                 /* We should free the atom name,
1288                  * but it might be used multiple times in the symtab.
1289                  * sfree(pdba->atomname[atind]);
1290                  */
1291                 for (int k = atind+1; k < pdba->nr; k++)
1292                 {
1293                     pdba->atom[k-1]     = pdba->atom[k];
1294                     pdba->atomname[k-1] = pdba->atomname[k];
1295                     copy_rvec(x[k], x[k-1]);
1296                 }
1297                 pdba->nr--;
1298                 bDeleted = true;
1299             }
1300         }
1301     }
1302
1303     return bDeleted;
1304 }
1305
1306 void match_atomnames_with_rtp(gmx::ArrayRef<PreprocessResidue>     usedPpResidues,
1307                               gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
1308                               t_atoms                             *pdba,
1309                               t_symtab                            *symtab,
1310                               gmx::ArrayRef<gmx::RVec>             x,
1311                               bool                                 bVerbose)
1312 {
1313     for (int i = 0; i < pdba->nr; i++)
1314     {
1315         const char        *oldnm           = *pdba->atomname[i];
1316         PreprocessResidue *localPpResidue  = &usedPpResidues[pdba->atom[i].resind];
1317         auto               found           = std::find_if(localPpResidue->atomname.begin(), localPpResidue->atomname.end(),
1318                                                           [&oldnm](char **name)
1319                                                           { return gmx::equalCaseInsensitive(oldnm, *name); });
1320         if (found == localPpResidue->atomname.end())
1321         {
1322             /* Not found yet, check if we have to rename this atom */
1323             if (match_atomnames_with_rtp_atom(pdba, x, symtab, i,
1324                                               localPpResidue, globalPatches[pdba->atom[i].resind],
1325                                               bVerbose))
1326             {
1327                 /* We deleted this atom, decrease the atom counter by 1. */
1328                 i--;
1329             }
1330         }
1331     }
1332 }
1333
1334 #define NUM_CMAP_ATOMS 5
1335 static void gen_cmap(InteractionTypeParameters *psb, gmx::ArrayRef<const PreprocessResidue> usedPpResidues, t_atoms *atoms)
1336 {
1337     int         residx;
1338     const char *ptr;
1339     t_resinfo  *resinfo = atoms->resinfo;
1340     int         nres    = atoms->nres;
1341     int         cmap_atomid[NUM_CMAP_ATOMS];
1342     int         cmap_chainnum = -1;
1343
1344     if (debug)
1345     {
1346         ptr = "cmap";
1347     }
1348     else
1349     {
1350         ptr = "check";
1351     }
1352
1353     fprintf(stderr, "Making cmap torsions...\n");
1354     int i = 0;
1355     /* Most cmap entries use the N atom from the next residue, so the last
1356      * residue should not have its CMAP entry in that case, but for things like
1357      * dipeptides we sometimes define a complete CMAP entry inside a residue,
1358      * and in this case we need to process everything through the last residue.
1359      */
1360     for (residx = 0; residx < nres; residx++)
1361     {
1362         /* Add CMAP terms from the list of CMAP interactions */
1363         for (const auto &b : usedPpResidues[residx].rb[ebtsCMAP].b)
1364         {
1365             bool bAddCMAP = true;
1366             /* Loop over atoms in a candidate CMAP interaction and
1367              * check that they exist, are from the same chain and are
1368              * from residues labelled as protein. */
1369             for (int k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1370             {
1371                 /* Assign the pointer to the name of the next reference atom.
1372                  * This can use -/+ labels to refer to previous/next residue.
1373                  */
1374                 const char *pname = b.a[k].c_str();
1375                 /* Skip this CMAP entry if it refers to residues before the
1376                  * first or after the last residue.
1377                  */
1378                 if (((strchr(pname, '-') != nullptr) && (residx == 0)) ||
1379                     ((strchr(pname, '+') != nullptr) && (residx == nres-1)))
1380                 {
1381                     bAddCMAP = false;
1382                     break;
1383                 }
1384
1385                 cmap_atomid[k] = search_atom(pname,
1386                                              i, atoms, ptr, TRUE);
1387                 bAddCMAP = bAddCMAP && (cmap_atomid[k] != -1);
1388                 if (!bAddCMAP)
1389                 {
1390                     /* This break is necessary, because cmap_atomid[k]
1391                      * == -1 cannot be safely used as an index
1392                      * into the atom array. */
1393                     break;
1394                 }
1395                 int this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1396                 if (0 == k)
1397                 {
1398                     cmap_chainnum = resinfo[this_residue_index].chainnum;
1399                 }
1400                 else
1401                 {
1402                     /* Does the residue for this atom have the same
1403                      * chain number as the residues for previous
1404                      * atoms? */
1405                     bAddCMAP = bAddCMAP &&
1406                         cmap_chainnum == resinfo[this_residue_index].chainnum;
1407                 }
1408                 /* Here we used to check that the residuetype was protein and
1409                  * disable bAddCMAP if that was not the case. However, some
1410                  * special residues (say, alanine dipeptides) might not adhere
1411                  * to standard naming, and if we start calling them normal
1412                  * protein residues the user will be bugged to select termini.
1413                  *
1414                  * Instead, I believe that the right course of action is to
1415                  * keep the CMAP interaction if it is present in the RTP file
1416                  * and we correctly identified all atoms (which is the case
1417                  * if we got here).
1418                  */
1419             }
1420
1421             if (bAddCMAP)
1422             {
1423                 add_cmap_param(psb, cmap_atomid[0], cmap_atomid[1], cmap_atomid[2], cmap_atomid[3], cmap_atomid[4], b.s.c_str());
1424             }
1425         }
1426
1427         if (residx < nres-1)
1428         {
1429             while (atoms->atom[i].resind < residx+1)
1430             {
1431                 i++;
1432             }
1433         }
1434     }
1435     /* Start the next residue */
1436 }
1437
1438 static void
1439 scrub_charge_groups(int *cgnr, int natoms)
1440 {
1441     int i;
1442
1443     for (i = 0; i < natoms; i++)
1444     {
1445         cgnr[i] = i+1;
1446     }
1447 }
1448
1449
1450 void pdb2top(FILE *top_file, const char *posre_fn, const char *molname,
1451              t_atoms *atoms,
1452              std::vector<gmx::RVec> *x, PreprocessingAtomTypes *atype, t_symtab *tab,
1453              gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
1454              gmx::ArrayRef<PreprocessResidue> usedPpResidues,
1455              gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
1456              bool bAllowMissing,
1457              bool bVsites, bool bVsiteAromatics,
1458              const char *ffdir,
1459              real mHmult,
1460              gmx::ArrayRef<const DisulfideBond> ssbonds,
1461              real long_bond_dist, real short_bond_dist,
1462              bool bDeuterate, bool bChargeGroups, bool bCmap,
1463              bool bRenumRes, bool bRTPresname)
1464 {
1465     std::array<InteractionTypeParameters, F_NRE> plist;
1466     t_excls                                     *excls;
1467     t_nextnb                                     nnb;
1468     int                                         *cgnr;
1469     int                                         *vsite_type;
1470     int                                          i, nmissat;
1471     int                                          bts[ebtsNR];
1472
1473     ResidueType rt;
1474
1475     /* Make bonds */
1476     at2bonds(&(plist[F_BONDS]), globalPatches,
1477              atoms, *x,
1478              long_bond_dist, short_bond_dist);
1479
1480     /* specbonds: disulphide bonds & heme-his */
1481     do_ssbonds(&(plist[F_BONDS]),
1482                atoms, ssbonds,
1483                bAllowMissing);
1484
1485     nmissat = name2type(atoms, &cgnr, usedPpResidues, &rt);
1486     if (nmissat)
1487     {
1488         if (bAllowMissing)
1489         {
1490             fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1491                     nmissat, molname);
1492         }
1493         else
1494         {
1495             gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1496                       nmissat, molname);
1497         }
1498     }
1499
1500     /* Cleanup bonds (sort and rm doubles) */
1501     clean_bonds(&(plist[F_BONDS]));
1502
1503     snew(vsite_type, atoms->nr);
1504     for (i = 0; i < atoms->nr; i++)
1505     {
1506         vsite_type[i] = NOTSET;
1507     }
1508     if (bVsites)
1509     {
1510         if (bVsiteAromatics)
1511         {
1512             fprintf(stdout, "The conversion of aromatic rings into virtual sites is deprecated "
1513                     "and may be removed in a future version of GROMACS");
1514         }
1515         /* determine which atoms will be vsites and add dummy masses
1516            also renumber atom numbers in plist[0..F_NRE]! */
1517         do_vsites(rtpFFDB, atype, atoms, tab, x, plist,
1518                   &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1519     }
1520
1521     /* Make Angles and Dihedrals */
1522     fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1523     snew(excls, atoms->nr);
1524     init_nnb(&nnb, atoms->nr, 4);
1525     gen_nnb(&nnb, plist);
1526     print_nnb(&nnb, "NNB");
1527     gen_pad(&nnb, atoms, usedPpResidues, plist, excls, globalPatches, bAllowMissing);
1528     done_nnb(&nnb);
1529
1530     /* Make CMAP */
1531     if (bCmap)
1532     {
1533         gen_cmap(&(plist[F_CMAP]), usedPpResidues, atoms);
1534         if (plist[F_CMAP].size() > 0)
1535         {
1536             fprintf(stderr, "There are %4zu cmap torsion pairs\n",
1537                     plist[F_CMAP].size());
1538         }
1539     }
1540
1541     /* set mass of all remaining hydrogen atoms */
1542     if (mHmult != 1.0)
1543     {
1544         do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1545     }
1546     sfree(vsite_type);
1547
1548     /* Cleanup bonds (sort and rm doubles) */
1549     /* clean_bonds(&(plist[F_BONDS]));*/
1550
1551     fprintf(stderr,
1552             "There are %4zu dihedrals, %4zu impropers, %4zu angles\n"
1553             "          %4zu pairs,     %4zu bonds and  %4zu virtual sites\n",
1554             plist[F_PDIHS].size(), plist[F_IDIHS].size(), plist[F_ANGLES].size(),
1555             plist[F_LJ14].size(), plist[F_BONDS].size(),
1556             plist[F_VSITE2].size() +
1557             plist[F_VSITE3].size() +
1558             plist[F_VSITE3FD].size() +
1559             plist[F_VSITE3FAD].size() +
1560             plist[F_VSITE3OUT].size() +
1561             plist[F_VSITE4FD].size() +
1562             plist[F_VSITE4FDN].size() );
1563
1564     print_sums(atoms, FALSE);
1565
1566     if (!bChargeGroups)
1567     {
1568         scrub_charge_groups(cgnr, atoms->nr);
1569     }
1570
1571     if (bRenumRes)
1572     {
1573         for (i = 0; i < atoms->nres; i++)
1574         {
1575             atoms->resinfo[i].nr = i + 1;
1576             atoms->resinfo[i].ic = ' ';
1577         }
1578     }
1579
1580     if (top_file)
1581     {
1582         fprintf(stderr, "Writing topology\n");
1583         /* We can copy the bonded types from the first restp,
1584          * since the types have to be identical for all residues in one molecule.
1585          */
1586         for (i = 0; i < ebtsNR; i++)
1587         {
1588             bts[i] = usedPpResidues[0].rb[i].type;
1589         }
1590         write_top(top_file, posre_fn, molname,
1591                   atoms, bRTPresname,
1592                   bts, plist, excls, atype, cgnr, usedPpResidues[0].nrexcl);
1593     }
1594
1595
1596     /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1597     sfree(cgnr);
1598     for (i = 0; i < atoms->nr; i++)
1599     {
1600         sfree(excls[i].e);
1601     }
1602     sfree(excls);
1603 }