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