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