src/gromacs/gmxpreprocess/pdb2gmx.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 "pdb2gmx.h"
  40
  41 #include <cctype>
  42 #include <cstdlib>
  43 #include <cstring>
  44 #include <ctime>
  45
  46 #include <algorithm>
  47 #include <memory>
  48 #include <string>
  49 #include <vector>
  50
  51 #include "gromacs/commandline/cmdlineoptionsmodule.h"
  52 #include "gromacs/fileio/confio.h"
  53 #include "gromacs/fileio/filetypes.h"
  54 #include "gromacs/fileio/gmxfio.h"
  55 #include "gromacs/fileio/pdbio.h"
  56 #include "gromacs/gmxlib/conformation-utilities.h"
  57 #include "gromacs/gmxpreprocess/fflibutil.h"
  58 #include "gromacs/gmxpreprocess/genhydro.h"
  59 #include "gromacs/gmxpreprocess/grompp-impl.h"
  60 #include "gromacs/gmxpreprocess/h_db.h"
  61 #include "gromacs/gmxpreprocess/hackblock.h"
  62 #include "gromacs/gmxpreprocess/hizzie.h"
  63 #include "gromacs/gmxpreprocess/pdb2top.h"
  64 #include "gromacs/gmxpreprocess/pgutil.h"
  65 #include "gromacs/gmxpreprocess/resall.h"
  66 #include "gromacs/gmxpreprocess/specbond.h"
  67 #include "gromacs/gmxpreprocess/ter_db.h"
  68 #include "gromacs/gmxpreprocess/toputil.h"
  69 #include "gromacs/gmxpreprocess/xlate.h"
  70 #include "gromacs/math/vec.h"
  71 #include "gromacs/options/basicoptions.h"
  72 #include "gromacs/options/filenameoption.h"
  73 #include "gromacs/options/ioptionscontainer.h"
  74 #include "gromacs/topology/atomprop.h"
  75 #include "gromacs/topology/block.h"
  76 #include "gromacs/topology/index.h"
  77 #include "gromacs/topology/residuetypes.h"
  78 #include "gromacs/topology/symtab.h"
  79 #include "gromacs/topology/topology.h"
  80 #include "gromacs/utility/dir_separator.h"
  81 #include "gromacs/utility/exceptions.h"
  82 #include "gromacs/utility/fatalerror.h"
  83 #include "gromacs/utility/path.h"
  84 #include "gromacs/utility/smalloc.h"
  85 #include "gromacs/utility/strdb.h"
  86 #include "gromacs/utility/stringutil.h"
  87
  88 #define RTP_MAXCHAR 5
  89 typedef struct {
  90     char gmx[RTP_MAXCHAR+2];
  91     char main[RTP_MAXCHAR+2];
  92     char nter[RTP_MAXCHAR+2];
  93     char cter[RTP_MAXCHAR+2];
  94     char bter[RTP_MAXCHAR+2];
  95 } rtprename_t;
  96
  97 static const char *res2bb_notermini(const char *name,
  98                                     int nrr, const rtprename_t *rr)
  99 {
 100     /* NOTE: This function returns the main building block name,
 101      *       it does not take terminal renaming into account.
 102      */
 103     int i;
 104
 105     i = 0;
 106     while (i < nrr && !gmx::equalCaseInsensitive(name, rr[i].gmx))
 107     {
 108         i++;
 109     }
 110
 111     return (i < nrr ? rr[i].main : name);
 112 }
 113
 114 static const char *select_res(int nr, int resnr,
 115                               const char *name[], const char *expl[],
 116                               const char *title,
 117                               int nrr, const rtprename_t *rr)
 118 {
 119     int sel = 0;
 120
 121     printf("Which %s type do you want for residue %d\n", title, resnr+1);
 122     for (sel = 0; (sel < nr); sel++)
 123     {
 124         printf("%d. %s (%s)\n",
 125                sel, expl[sel], res2bb_notermini(name[sel], nrr, rr));
 126     }
 127     printf("\nType a number:"); fflush(stdout);
 128
 129     if (scanf("%d", &sel) != 1)
 130     {
 131         gmx_fatal(FARGS, "Answer me for res %s %d!", title, resnr+1);
 132     }
 133
 134     return name[sel];
 135 }
 136
 137 static const char *get_asptp(int resnr, int nrr, const rtprename_t *rr)
 138 {
 139     enum {
 140         easp, easpH, easpNR
 141     };
 142     const char *lh[easpNR]   = { "ASP", "ASPH" };
 143     const char *expl[easpNR] = {
 144         "Not protonated (charge -1)",
 145         "Protonated (charge 0)"
 146     };
 147
 148     return select_res(easpNR, resnr, lh, expl, "ASPARTIC ACID", nrr, rr);
 149 }
 150
 151 static const char *get_glutp(int resnr, int nrr, const rtprename_t *rr)
 152 {
 153     enum {
 154         eglu, egluH, egluNR
 155     };
 156     const char *lh[egluNR]   = { "GLU", "GLUH" };
 157     const char *expl[egluNR] = {
 158         "Not protonated (charge -1)",
 159         "Protonated (charge 0)"
 160     };
 161
 162     return select_res(egluNR, resnr, lh, expl, "GLUTAMIC ACID", nrr, rr);
 163 }
 164
 165 static const char *get_glntp(int resnr, int nrr, const rtprename_t *rr)
 166 {
 167     enum {
 168         egln, eglnH, eglnNR
 169     };
 170     const char *lh[eglnNR]   = { "GLN", "QLN" };
 171     const char *expl[eglnNR] = {
 172         "Not protonated (charge 0)",
 173         "Protonated (charge +1)"
 174     };
 175
 176     return select_res(eglnNR, resnr, lh, expl, "GLUTAMINE", nrr, rr);
 177 }
 178
 179 static const char *get_lystp(int resnr, int nrr, const rtprename_t *rr)
 180 {
 181     enum {
 182         elys, elysH, elysNR
 183     };
 184     const  char *lh[elysNR]   = { "LYSN", "LYS" };
 185     const char  *expl[elysNR] = {
 186         "Not protonated (charge 0)",
 187         "Protonated (charge +1)"
 188     };
 189
 190     return select_res(elysNR, resnr, lh, expl, "LYSINE", nrr, rr);
 191 }
 192
 193 static const char *get_argtp(int resnr, int nrr, const rtprename_t *rr)
 194 {
 195     enum {
 196         earg, eargH, eargNR
 197     };
 198     const  char *lh[eargNR]   = { "ARGN", "ARG" };
 199     const char  *expl[eargNR] = {
 200         "Not protonated (charge 0)",
 201         "Protonated (charge +1)"
 202     };
 203
 204     return select_res(eargNR, resnr, lh, expl, "ARGININE", nrr, rr);
 205 }
 206
 207 static const char *get_histp(int resnr, int nrr, const rtprename_t *rr)
 208 {
 209     const char *expl[ehisNR] = {
 210         "H on ND1 only",
 211         "H on NE2 only",
 212         "H on ND1 and NE2",
 213         "Coupled to Heme"
 214     };
 215
 216     return select_res(ehisNR, resnr, hh, expl, "HISTIDINE", nrr, rr);
 217 }
 218
 219 static void read_rtprename(const char *fname, FILE *fp,
 220                            int *nrtprename, rtprename_t **rtprename)
 221 {
 222     char         line[STRLEN], buf[STRLEN];
 223     int          n;
 224     rtprename_t *rr;
 225     int          ncol, nc;
 226
 227     n  = *nrtprename;
 228     rr = *rtprename;
 229
 230     ncol = 0;
 231     while (get_a_line(fp, line, STRLEN))
 232     {
 233         srenew(rr, n+1);
 234         /* line is NULL-terminated and length<STRLEN, so final arg cannot overflow.
 235          * For other args, we read up to 6 chars (so we can detect if the length is > 5).
 236          * Note that the buffer length has been increased to 7 to allow this,
 237          * so we just need to make sure the strings have zero-length initially.
 238          */
 239         rr[n].gmx[0]  = '\0';
 240         rr[n].main[0] = '\0';
 241         rr[n].nter[0] = '\0';
 242         rr[n].cter[0] = '\0';
 243         rr[n].bter[0] = '\0';
 244         nc            = sscanf(line, "%6s %6s %6s %6s %6s %s",
 245                                rr[n].gmx, rr[n].main, rr[n].nter, rr[n].cter, rr[n].bter, buf);
 246         if (strlen(rr[n].gmx) > RTP_MAXCHAR || strlen(rr[n].main) > RTP_MAXCHAR ||
 247             strlen(rr[n].nter) > RTP_MAXCHAR || strlen(rr[n].cter) > RTP_MAXCHAR || strlen(rr[n].bter) > RTP_MAXCHAR)
 248         {
 249             gmx_fatal(FARGS, "Residue renaming database '%s' has strings longer than %d chars in first 5 columns:\n%s",
 250                       fname, RTP_MAXCHAR, line);
 251         }
 252
 253         if (ncol == 0)
 254         {
 255             if (nc != 2 && nc != 5)
 256             {
 257                 gmx_fatal(FARGS, "Residue renaming database '%s' has %d columns instead of %d or %d", fname, ncol, 2, 5);
 258             }
 259             ncol = nc;
 260         }
 261         else if (nc != ncol)
 262         {
 263             gmx_fatal(FARGS, "A line in residue renaming database '%s' has %d columns, while previous lines have %d columns", fname, nc, ncol);
 264         }
 265
 266         if (nc == 2)
 267         {
 268             /* This file does not have special termini names, copy them from main */
 269             strcpy(rr[n].nter, rr[n].main);
 270             strcpy(rr[n].cter, rr[n].main);
 271             strcpy(rr[n].bter, rr[n].main);
 272         }
 273
 274         n++;
 275     }
 276
 277     *nrtprename = n;
 278     *rtprename  = rr;
 279 }
 280
 281 static char *search_resrename(int nrr, rtprename_t *rr,
 282                               const char *name,
 283                               bool bStart, bool bEnd,
 284                               bool bCompareFFRTPname)
 285 {
 286     char *nn;
 287     int   i;
 288
 289     nn = nullptr;
 290
 291     i = 0;
 292     while (i < nrr && ((!bCompareFFRTPname && strcmp(name, rr[i].gmx)  != 0) ||
 293                        ( bCompareFFRTPname && strcmp(name, rr[i].main) != 0)))
 294     {
 295         i++;
 296     }
 297
 298     /* If found in the database, rename this residue's rtp building block,
 299      * otherwise keep the old name.
 300      */
 301     if (i < nrr)
 302     {
 303         if (bStart && bEnd)
 304         {
 305             nn = rr[i].bter;
 306         }
 307         else if (bStart)
 308         {
 309             nn = rr[i].nter;
 310         }
 311         else if (bEnd)
 312         {
 313             nn = rr[i].cter;
 314         }
 315         else
 316         {
 317             nn = rr[i].main;
 318         }
 319
 320         if (nn[0] == '-')
 321         {
 322             gmx_fatal(FARGS, "In the chosen force field there is no residue type for '%s'%s", name, bStart ? ( bEnd ? " as a standalone (starting & ending) residue" : " as a starting terminus") : (bEnd ? " as an ending terminus" : ""));
 323         }
 324     }
 325
 326     return nn;
 327 }
 328
 329 static void rename_resrtp(t_atoms *pdba, int nterpairs, const int *r_start, const int *r_end,
 330                           int nrr, rtprename_t *rr, t_symtab *symtab,
 331                           bool bVerbose)
 332 {
 333     int      r, j;
 334     bool     bStart, bEnd;
 335     char    *nn;
 336     bool     bFFRTPTERRNM;
 337
 338     bFFRTPTERRNM = (getenv("GMX_NO_FFRTP_TER_RENAME") == nullptr);
 339
 340     for (r = 0; r < pdba->nres; r++)
 341     {
 342         bStart = false;
 343         bEnd   = false;
 344         for (j = 0; j < nterpairs; j++)
 345         {
 346             if (r == r_start[j])
 347             {
 348                 bStart = true;
 349             }
 350         }
 351         for (j = 0; j < nterpairs; j++)
 352         {
 353             if (r == r_end[j])
 354             {
 355                 bEnd = true;
 356             }
 357         }
 358
 359         nn = search_resrename(nrr, rr, *pdba->resinfo[r].rtp, bStart, bEnd, false);
 360
 361         if (bFFRTPTERRNM && nn == nullptr && (bStart || bEnd))
 362         {
 363             /* This is a terminal residue, but the residue name,
 364              * currently stored in .rtp, is not a standard residue name,
 365              * but probably a force field specific rtp name.
 366              * Check if we need to rename it because it is terminal.
 367              */
 368             nn = search_resrename(nrr, rr,
 369                                   *pdba->resinfo[r].rtp, bStart, bEnd, true);
 370         }
 371
 372         if (nn != nullptr && strcmp(*pdba->resinfo[r].rtp, nn) != 0)
 373         {
 374             if (bVerbose)
 375             {
 376                 printf("Changing rtp entry of residue %d %s to '%s'\n",
 377                        pdba->resinfo[r].nr, *pdba->resinfo[r].name, nn);
 378             }
 379             pdba->resinfo[r].rtp = put_symtab(symtab, nn);
 380         }
 381     }
 382 }
 383
 384 static void pdbres_to_gmxrtp(t_atoms *pdba)
 385 {
 386     int i;
 387
 388     for (i = 0; (i < pdba->nres); i++)
 389     {
 390         if (pdba->resinfo[i].rtp == nullptr)
 391         {
 392             pdba->resinfo[i].rtp = pdba->resinfo[i].name;
 393         }
 394     }
 395 }
 396
 397 static void rename_pdbres(t_atoms *pdba, const char *oldnm, const char *newnm,
 398                           bool bFullCompare, t_symtab *symtab)
 399 {
 400     char *resnm;
 401     int   i;
 402
 403     for (i = 0; (i < pdba->nres); i++)
 404     {
 405         resnm = *pdba->resinfo[i].name;
 406         if ((bFullCompare && (gmx::equalCaseInsensitive(resnm, oldnm))) ||
 407             (!bFullCompare && strstr(resnm, oldnm) != nullptr))
 408         {
 409             /* Rename the residue name (not the rtp name) */
 410             pdba->resinfo[i].name = put_symtab(symtab, newnm);
 411         }
 412     }
 413 }
 414
 415 static void rename_bb(t_atoms *pdba, const char *oldnm, const char *newnm,
 416                       bool bFullCompare, t_symtab *symtab)
 417 {
 418     char *bbnm;
 419     int   i;
 420
 421     for (i = 0; (i < pdba->nres); i++)
 422     {
 423         /* We have not set the rtp name yes, use the residue name */
 424         bbnm = *pdba->resinfo[i].name;
 425         if ((bFullCompare && (gmx::equalCaseInsensitive(bbnm, oldnm))) ||
 426             (!bFullCompare && strstr(bbnm, oldnm) != nullptr))
 427         {
 428             /* Change the rtp builing block name */
 429             pdba->resinfo[i].rtp = put_symtab(symtab, newnm);
 430         }
 431     }
 432 }
 433
 434 static void rename_bbint(t_atoms *pdba, const char *oldnm,
 435                          const char *gettp(int, int, const rtprename_t *),
 436                          bool bFullCompare,
 437                          t_symtab *symtab,
 438                          int nrr, const rtprename_t *rr)
 439 {
 440     int         i;
 441     const char *ptr;
 442     char       *bbnm;
 443
 444     for (i = 0; i < pdba->nres; i++)
 445     {
 446         /* We have not set the rtp name yet, use the residue name */
 447         bbnm = *pdba->resinfo[i].name;
 448         if ((bFullCompare && (strcmp(bbnm, oldnm) == 0)) ||
 449             (!bFullCompare && strstr(bbnm, oldnm) != nullptr))
 450         {
 451             ptr                  = gettp(i, nrr, rr);
 452             pdba->resinfo[i].rtp = put_symtab(symtab, ptr);
 453         }
 454     }
 455 }
 456
 457 static void check_occupancy(t_atoms *atoms, const char *filename, bool bVerbose)
 458 {
 459     int i, ftp;
 460     int nzero   = 0;
 461     int nnotone = 0;
 462
 463     ftp = fn2ftp(filename);
 464     if (!atoms->pdbinfo || ((ftp != efPDB) && (ftp != efBRK) && (ftp != efENT)))
 465     {
 466         fprintf(stderr, "No occupancies in %s\n", filename);
 467     }
 468     else
 469     {
 470         for (i = 0; (i < atoms->nr); i++)
 471         {
 472             if (atoms->pdbinfo[i].occup != 1)
 473             {
 474                 if (bVerbose)
 475                 {
 476                     fprintf(stderr, "Occupancy for atom %s%d-%s is %f rather than 1\n",
 477                             *atoms->resinfo[atoms->atom[i].resind].name,
 478                             atoms->resinfo[ atoms->atom[i].resind].nr,
 479                             *atoms->atomname[i],
 480                             atoms->pdbinfo[i].occup);
 481                 }
 482                 if (atoms->pdbinfo[i].occup == 0)
 483                 {
 484                     nzero++;
 485                 }
 486                 else
 487                 {
 488                     nnotone++;
 489                 }
 490             }
 491         }
 492         if (nzero == atoms->nr)
 493         {
 494             fprintf(stderr, "All occupancy fields zero. This is probably not an X-Ray structure\n");
 495         }
 496         else if ((nzero > 0) || (nnotone > 0))
 497         {
 498             fprintf(stderr,
 499                     "\n"
 500                     "WARNING: there were %d atoms with zero occupancy and %d atoms with\n"
 501                     "         occupancy unequal to one (out of %d atoms). Check your pdb file.\n"
 502                     "\n",
 503                     nzero, nnotone, atoms->nr);
 504         }
 505         else
 506         {
 507             fprintf(stderr, "All occupancies are one\n");
 508         }
 509     }
 510 }
 511
 512 static void write_posres(const char *fn, t_atoms *pdba, real fc)
 513 {
 514     FILE *fp;
 515     int   i;
 516
 517     fp = gmx_fio_fopen(fn, "w");
 518     fprintf(fp,
 519             "; In this topology include file, you will find position restraint\n"
 520             "; entries for all the heavy atoms in your original pdb file.\n"
 521             "; This means that all the protons which were added by pdb2gmx are\n"
 522             "; not restrained.\n"
 523             "\n"
 524             "[ position_restraints ]\n"
 525             "; %4s%6s%8s%8s%8s\n", "atom", "type", "fx", "fy", "fz"
 526             );
 527     for (i = 0; (i < pdba->nr); i++)
 528     {
 529         if (!is_hydrogen(*pdba->atomname[i]) && !is_dummymass(*pdba->atomname[i]))
 530         {
 531             fprintf(fp, "%6d%6d  %g  %g  %g\n", i+1, 1, fc, fc, fc);
 532         }
 533     }
 534     gmx_fio_fclose(fp);
 535 }
 536
 537 static int read_pdball(const char *inf, bool bOutput, const char *outf, char **title,
 538                        t_atoms *atoms, rvec **x,
 539                        int *ePBC, matrix box, bool bRemoveH,
 540                        t_symtab *symtab, ResidueType *rt, const char *watres,
 541                        AtomProperties *aps, bool bVerbose)
 542 /* Read a pdb file. (containing proteins) */
 543 {
 544     int natom, new_natom, i;
 545
 546     /* READ IT */
 547     printf("Reading %s...\n", inf);
 548     readConfAndAtoms(inf, symtab, title, atoms, ePBC, x, nullptr, box);
 549     natom           = atoms->nr;
 550     if (atoms->pdbinfo == nullptr)
 551     {
 552         snew(atoms->pdbinfo, atoms->nr);
 553     }
 554     if (fn2ftp(inf) == efPDB)
 555     {
 556         get_pdb_atomnumber(atoms, aps);
 557     }
 558     if (bRemoveH)
 559     {
 560         new_natom = 0;
 561         for (i = 0; i < atoms->nr; i++)
 562         {
 563             if (!is_hydrogen(*atoms->atomname[i]))
 564             {
 565                 atoms->atom[new_natom]     = atoms->atom[i];
 566                 atoms->atomname[new_natom] = atoms->atomname[i];
 567                 atoms->pdbinfo[new_natom]  = atoms->pdbinfo[i];
 568                 copy_rvec((*x)[i], (*x)[new_natom]);
 569                 new_natom++;
 570             }
 571         }
 572         atoms->nr = new_natom;
 573         natom     = new_natom;
 574     }
 575
 576     printf("Read");
 577     if (title[0])
 578     {
 579         printf(" '%s',", *title);
 580     }
 581     printf(" %d atoms\n", natom);
 582
 583     /* Rename residues */
 584     rename_pdbres(atoms, "HOH", watres, false, symtab);
 585     rename_pdbres(atoms, "SOL", watres, false, symtab);
 586     rename_pdbres(atoms, "WAT", watres, false, symtab);
 587
 588     rename_atoms("xlateat.dat", nullptr,
 589                  atoms, symtab, nullptr, true,
 590                  rt, true, bVerbose);
 591
 592     if (natom == 0)
 593     {
 594         return 0;
 595     }
 596     if (bOutput)
 597     {
 598         write_sto_conf(outf, *title, atoms, *x, nullptr, *ePBC, box);
 599     }
 600
 601     return natom;
 602 }
 603
 604 static void process_chain(t_atoms *pdba, rvec *x,
 605                           bool bTrpU, bool bPheU, bool bTyrU,
 606                           bool bLysMan, bool bAspMan, bool bGluMan,
 607                           bool bHisMan, bool bArgMan, bool bGlnMan,
 608                           real angle, real distance, t_symtab *symtab,
 609                           int nrr, const rtprename_t *rr)
 610 {
 611     /* Rename aromatics, lys, asp and histidine */
 612     if (bTyrU)
 613     {
 614         rename_bb(pdba, "TYR", "TYRU", false, symtab);
 615     }
 616     if (bTrpU)
 617     {
 618         rename_bb(pdba, "TRP", "TRPU", false, symtab);
 619     }
 620     if (bPheU)
 621     {
 622         rename_bb(pdba, "PHE", "PHEU", false, symtab);
 623     }
 624     if (bLysMan)
 625     {
 626         rename_bbint(pdba, "LYS", get_lystp, false, symtab, nrr, rr);
 627     }
 628     if (bArgMan)
 629     {
 630         rename_bbint(pdba, "ARG", get_argtp, false, symtab, nrr, rr);
 631     }
 632     if (bGlnMan)
 633     {
 634         rename_bbint(pdba, "GLN", get_glntp, false, symtab, nrr, rr);
 635     }
 636     if (bAspMan)
 637     {
 638         rename_bbint(pdba, "ASP", get_asptp, false, symtab, nrr, rr);
 639     }
 640     else
 641     {
 642         rename_bb(pdba, "ASPH", "ASP", false, symtab);
 643     }
 644     if (bGluMan)
 645     {
 646         rename_bbint(pdba, "GLU", get_glutp, false, symtab, nrr, rr);
 647     }
 648     else
 649     {
 650         rename_bb(pdba, "GLUH", "GLU", false, symtab);
 651     }
 652
 653     if (!bHisMan)
 654     {
 655         set_histp(pdba, x, angle, distance);
 656     }
 657     else
 658     {
 659         rename_bbint(pdba, "HIS", get_histp, true, symtab, nrr, rr);
 660     }
 661
 662     /* Initialize the rtp builing block names with the residue names
 663      * for the residues that have not been processed above.
 664      */
 665     pdbres_to_gmxrtp(pdba);
 666
 667     /* Now we have all rtp names set.
 668      * The rtp names will conform to Gromacs naming,
 669      * unless the input pdb file contained one or more force field specific
 670      * rtp names as residue names.
 671      */
 672 }
 673
 674 /* struct for sorting the atoms from the pdb file */
 675 typedef struct {
 676     int  resnr;  /* residue number               */
 677     int  j;      /* database order index         */
 678     int  index;  /* original atom number         */
 679     char anm1;   /* second letter of atom name   */
 680     char altloc; /* alternate location indicator */
 681 } t_pdbindex;
 682
 683 static bool pdbicomp(const t_pdbindex &a, const t_pdbindex &b)
 684 {
 685     int d = (a.resnr - b.resnr);
 686     if (d == 0)
 687     {
 688         d = (a.j - b.j);
 689         if (d == 0)
 690         {
 691             d = (a.anm1 - b.anm1);
 692             if (d == 0)
 693             {
 694                 d = (a.altloc - b.altloc);
 695             }
 696         }
 697     }
 698     return d < 0;
 699 }
 700
 701 static void sort_pdbatoms(t_restp restp[],
 702                           int natoms, t_atoms **pdbaptr, rvec **x,
 703                           t_blocka *block, char ***gnames)
 704 {
 705     t_atoms     *pdba, *pdbnew;
 706     rvec       **xnew;
 707     t_restp     *rptr;
 708     t_pdbindex  *pdbi;
 709     char        *atomnm;
 710
 711     pdba   = *pdbaptr;
 712     natoms = pdba->nr;
 713     pdbnew = nullptr;
 714     snew(xnew, 1);
 715     snew(pdbi, natoms);
 716
 717     for (int i = 0; i < natoms; i++)
 718     {
 719         atomnm = *pdba->atomname[i];
 720         rptr   = &restp[pdba->atom[i].resind];
 721         int j = std::find_if(rptr->atomname, rptr->atomname+rptr->natom,
 722                              [&atomnm](char** it){return gmx::equalCaseInsensitive(atomnm, *it); })
 723             - rptr->atomname;
 724         if (j == rptr->natom)
 725         {
 726             char buf[STRLEN];
 727
 728             sprintf(buf,
 729                     "Atom %s in residue %s %d was not found in rtp entry %s with %d atoms\n"
 730                     "while sorting atoms.\n%s", atomnm,
 731                     *pdba->resinfo[pdba->atom[i].resind].name,
 732                     pdba->resinfo[pdba->atom[i].resind].nr,
 733                     rptr->resname,
 734                     rptr->natom,
 735                     is_hydrogen(atomnm) ?
 736                     "\nFor a hydrogen, this can be a different protonation state, or it\n"
 737                     "might have had a different number in the PDB file and was rebuilt\n"
 738                     "(it might for instance have been H3, and we only expected H1 & H2).\n"
 739                     "Note that hydrogens might have been added to the entry for the N-terminus.\n"
 740                     "Remove this hydrogen or choose a different protonation state to solve it.\n"
 741                     "Option -ignh will ignore all hydrogens in the input." : ".");
 742             gmx_fatal(FARGS, "%s", buf);
 743         }
 744         /* make shadow array to be sorted into indexgroup */
 745         pdbi[i].resnr  = pdba->atom[i].resind;
 746         pdbi[i].j      = j;
 747         pdbi[i].index  = i;
 748         pdbi[i].anm1   = atomnm[1];
 749         pdbi[i].altloc = pdba->pdbinfo[i].altloc;
 750     }
 751     std::sort(pdbi, pdbi+natoms, pdbicomp);
 752
 753     /* pdba is sorted in pdbnew using the pdbi index */
 754     std::vector<int> a(natoms);
 755     snew(pdbnew, 1);
 756     init_t_atoms(pdbnew, natoms, true);
 757     snew(*xnew, natoms);
 758     pdbnew->nr   = pdba->nr;
 759     pdbnew->nres = pdba->nres;
 760     sfree(pdbnew->resinfo);
 761     pdbnew->resinfo = pdba->resinfo;
 762     for (int i = 0; i < natoms; i++)
 763     {
 764         pdbnew->atom[i]     = pdba->atom[pdbi[i].index];
 765         pdbnew->atomname[i] = pdba->atomname[pdbi[i].index];
 766         pdbnew->pdbinfo[i]  = pdba->pdbinfo[pdbi[i].index];
 767         copy_rvec((*x)[pdbi[i].index], (*xnew)[i]);
 768         /* make indexgroup in block */
 769         a[i] = pdbi[i].index;
 770     }
 771     /* clean up */
 772     sfree(pdba->atomname);
 773     sfree(pdba->atom);
 774     sfree(pdba->pdbinfo);
 775     sfree(pdba);
 776     sfree(*x);
 777     /* copy the sorted pdbnew back to pdba */
 778     *pdbaptr = pdbnew;
 779     *x       = *xnew;
 780     add_grp(block, gnames, a, "prot_sort");
 781     sfree(xnew);
 782     sfree(pdbi);
 783 }
 784
 785 static int remove_duplicate_atoms(t_atoms *pdba, rvec x[], bool bVerbose)
 786 {
 787     int        i, j, oldnatoms, ndel;
 788     t_resinfo *ri;
 789
 790     printf("Checking for duplicate atoms....\n");
 791     oldnatoms    = pdba->nr;
 792     ndel         = 0;
 793     /* NOTE: pdba->nr is modified inside the loop */
 794     for (i = 1; (i < pdba->nr); i++)
 795     {
 796         /* compare 'i' and 'i-1', throw away 'i' if they are identical
 797            this is a 'while' because multiple alternate locations can be present */
 798         while ( (i < pdba->nr) &&
 799                 (pdba->atom[i-1].resind == pdba->atom[i].resind) &&
 800                 (strcmp(*pdba->atomname[i-1], *pdba->atomname[i]) == 0) )
 801         {
 802             ndel++;
 803             if (bVerbose)
 804             {
 805                 ri = &pdba->resinfo[pdba->atom[i].resind];
 806                 printf("deleting duplicate atom %4s  %s%4d%c",
 807                        *pdba->atomname[i], *ri->name, ri->nr, ri->ic);
 808                 if (ri->chainid && (ri->chainid != ' '))
 809                 {
 810                     printf(" ch %c", ri->chainid);
 811                 }
 812                 if (pdba->pdbinfo)
 813                 {
 814                     if (pdba->pdbinfo[i].atomnr)
 815                     {
 816                         printf("  pdb nr %4d", pdba->pdbinfo[i].atomnr);
 817                     }
 818                     if (pdba->pdbinfo[i].altloc && (pdba->pdbinfo[i].altloc != ' '))
 819                     {
 820                         printf("  altloc %c", pdba->pdbinfo[i].altloc);
 821                     }
 822                 }
 823                 printf("\n");
 824             }
 825             pdba->nr--;
 826             /* We can not free, since it might be in the symtab */
 827             /* sfree(pdba->atomname[i]); */
 828             for (j = i; j < pdba->nr; j++)
 829             {
 830                 pdba->atom[j]     = pdba->atom[j+1];
 831                 pdba->atomname[j] = pdba->atomname[j+1];
 832                 if (pdba->pdbinfo)
 833                 {
 834                     pdba->pdbinfo[j]  = pdba->pdbinfo[j+1];
 835                 }
 836                 copy_rvec(x[j+1], x[j]);
 837             }
 838             srenew(pdba->atom,     pdba->nr);
 839             /* srenew(pdba->atomname, pdba->nr); */
 840             srenew(pdba->pdbinfo,  pdba->nr);
 841         }
 842     }
 843     if (pdba->nr != oldnatoms)
 844     {
 845         printf("Now there are %d atoms. Deleted %d duplicates.\n", pdba->nr, ndel);
 846     }
 847
 848     return pdba->nr;
 849 }
 850
 851 static void
 852 checkResidueTypeSanity(t_atoms     *pdba,
 853                        int          r0,
 854                        int          r1,
 855                        ResidueType *rt)
 856 {
 857     std::string startResidueString = gmx::formatString("%s%d", *pdba->resinfo[r0].name, pdba->resinfo[r0].nr);
 858     std::string endResidueString   = gmx::formatString("%s%d", *pdba->resinfo[r1-1].name, pdba->resinfo[r1-1].nr);
 859
 860     // Check whether all residues in chain have the same chain ID.
 861     bool         allResiduesHaveSameChainID = true;
 862     char         chainID0                   = pdba->resinfo[r0].chainid;
 863     char         chainID;
 864     std::string  residueString;
 865
 866     for (int i = r0 + 1; i < r1; i++)
 867     {
 868         chainID = pdba->resinfo[i].chainid;
 869         if (chainID != chainID0)
 870         {
 871             allResiduesHaveSameChainID  = false;
 872             residueString               = gmx::formatString("%s%d", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
 873             break;
 874         }
 875     }
 876
 877     if (!allResiduesHaveSameChainID)
 878     {
 879         gmx_fatal(FARGS,
 880                   "The chain covering the range %s--%s does not have a consistent chain ID. "
 881                   "The first residue has ID '%c', while residue %s has ID '%c'.",
 882                   startResidueString.c_str(), endResidueString.c_str(),
 883                   chainID0, residueString.c_str(), chainID);
 884     }
 885
 886     // At this point all residues have the same ID. If they are also non-blank
 887     // we can be a bit more aggressive and require the types match too.
 888     if (chainID0 != ' ')
 889     {
 890         bool        allResiduesHaveSameType = true;
 891         std::string restype;
 892         std::string restype0 = rt->typeNameForIndexedResidue(*pdba->resinfo[r0].name);
 893
 894         for (int i = r0 + 1; i < r1; i++)
 895         {
 896             restype = rt->typeNameForIndexedResidue(*pdba->resinfo[i].name);
 897             if (!gmx::equalCaseInsensitive(restype, restype0))
 898             {
 899                 allResiduesHaveSameType = false;
 900                 residueString           = gmx::formatString("%s%d", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
 901                 break;
 902             }
 903         }
 904
 905         if (!allResiduesHaveSameType)
 906         {
 907             gmx_fatal(FARGS,
 908                       "The residues in the chain %s--%s do not have a consistent type. "
 909                       "The first residue has type '%s', while residue %s is of type '%s'. "
 910                       "Either there is a mistake in your chain, or it includes nonstandard "
 911                       "residue names that have not yet been added to the residuetypes.dat "
 912                       "file in the GROMACS library directory. If there are other molecules "
 913                       "such as ligands, they should not have the same chain ID as the "
 914                       "adjacent protein chain since it's a separate molecule.",
 915                       startResidueString.c_str(), endResidueString.c_str(),
 916                       restype0.c_str(), residueString.c_str(), restype.c_str());
 917         }
 918     }
 919 }
 920
 921 static void find_nc_ter(t_atoms *pdba, int r0, int r1, int *r_start, int *r_end,
 922                         ResidueType *rt)
 923 {
 924     int         i;
 925     std::string p_startrestype;
 926
 927     *r_start = -1;
 928     *r_end   = -1;
 929
 930     int startWarnings = 0;
 931     int endWarnings   = 0;
 932     int ionNotes      = 0;
 933
 934     // Check that all residues have the same chain identifier, and if it is
 935     // non-blank we also require the residue types to match.
 936     checkResidueTypeSanity(pdba, r0, r1, rt);
 937
 938     // If we return correctly from checkResidueTypeSanity(), the only
 939     // remaining cases where we can have non-matching residue types is if
 940     // the chain ID was blank, which could be the case e.g. for a structure
 941     // read from a GRO file or other file types without chain information.
 942     // In that case we need to be a bit more liberal and detect chains based
 943     // on the residue type.
 944
 945     // If we get here, the chain ID must be identical for all residues
 946     char chainID = pdba->resinfo[r0].chainid;
 947
 948     /* Find the starting terminus (typially N or 5') */
 949     for (i = r0; i < r1 && *r_start == -1; i++)
 950     {
 951         p_startrestype = rt->typeNameForIndexedResidue(*pdba->resinfo[i].name);
 952         if (gmx::equalCaseInsensitive(p_startrestype, "Protein") ||
 953             gmx::equalCaseInsensitive(p_startrestype, "DNA") ||
 954             gmx::equalCaseInsensitive(p_startrestype, "RNA") )
 955         {
 956             printf("Identified residue %s%d as a starting terminus.\n", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
 957             *r_start = i;
 958         }
 959         else if (gmx::equalCaseInsensitive(p_startrestype, "Ion"))
 960         {
 961             if (ionNotes < 5)
 962             {
 963                 printf("Residue %s%d has type 'Ion', assuming it is not linked into a chain.\n", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
 964             }
 965             if (ionNotes == 4)
 966             {
 967                 printf("Disabling further notes about ions.\n");
 968             }
 969             ionNotes++;
 970         }
 971         else
 972         {
 973             if (startWarnings < 5)
 974             {
 975                 if (chainID == ' ')
 976                 {
 977                     printf("\nWarning: Starting residue %s%d in chain not identified as Protein/RNA/DNA.\n"
 978                            "This chain lacks identifiers, which makes it impossible to do strict\n"
 979                            "classification of the start/end residues. Here we need to guess this residue\n"
 980                            "should not be part of the chain and instead introduce a break, but that will\n"
 981                            "be catastrophic if they should in fact be linked. Please check your structure,\n"
 982                            "and add %s to residuetypes.dat if this was not correct.\n\n",
 983                            *pdba->resinfo[i].name, pdba->resinfo[i].nr, *pdba->resinfo[i].name);
 984                 }
 985                 else
 986                 {
 987                     printf("\nWarning: No residues in chain starting at %s%d identified as Protein/RNA/DNA.\n"
 988                            "This makes it impossible to link them into a molecule, which could either be\n"
 989                            "correct or a catastrophic error. Please check your structure, and add all\n"
 990                            "necessary residue names to residuetypes.dat if this was not correct.\n\n",
 991                            *pdba->resinfo[i].name, pdba->resinfo[i].nr);
 992                 }
 993             }
 994             if (startWarnings == 4)
 995             {
 996                 printf("Disabling further warnings about unidentified residues at start of chain.\n");
 997             }
 998             startWarnings++;
 999         }
1000     }
1001
1002     if (*r_start >= 0)
1003     {
1004         /* Go through the rest of the residues, check that they are the same class, and identify the ending terminus. */
1005         for (int i = *r_start; i < r1; i++)
1006         {
1007             std::string p_restype = rt->typeNameForIndexedResidue(*pdba->resinfo[i].name);
1008             if (gmx::equalCaseInsensitive(p_restype, p_startrestype) && endWarnings == 0)
1009             {
1010                 *r_end = i;
1011             }
1012             else if (gmx::equalCaseInsensitive(p_startrestype, "Ion"))
1013             {
1014                 if (ionNotes < 5)
1015                 {
1016                     printf("Residue %s%d has type 'Ion', assuming it is not linked into a chain.\n", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
1017                 }
1018                 if (ionNotes == 4)
1019                 {
1020                     printf("Disabling further notes about ions.\n");
1021                 }
1022                 ionNotes++;
1023             }
1024             else
1025             {
1026                 // This can only trigger if the chain ID is blank - otherwise the
1027                 // call to checkResidueTypeSanity() will have caught the problem.
1028                 if (endWarnings < 5)
1029                 {
1030                     printf("\nWarning: Residue %s%d in chain has different type ('%s') from\n"
1031                            "residue %s%d ('%s'). This chain lacks identifiers, which makes\n"
1032                            "it impossible to do strict classification of the start/end residues. Here we\n"
1033                            "need to guess this residue should not be part of the chain and instead\n"
1034                            "introduce a break, but that will be catastrophic if they should in fact be\n"
1035                            "linked. Please check your structure, and add %s to residuetypes.dat\n"
1036                            "if this was not correct.\n\n",
1037                            *pdba->resinfo[i].name, pdba->resinfo[i].nr, p_restype.c_str(),
1038                            *pdba->resinfo[*r_start].name, pdba->resinfo[*r_start].nr, p_startrestype.c_str(), *pdba->resinfo[i].name);
1039                 }
1040                 if (endWarnings == 4)
1041                 {
1042                     printf("Disabling further warnings about unidentified residues at end of chain.\n");
1043                 }
1044                 endWarnings++;
1045             }
1046         }
1047     }
1048
1049     if (*r_end >= 0)
1050     {
1051         printf("Identified residue %s%d as a ending terminus.\n", *pdba->resinfo[*r_end].name, pdba->resinfo[*r_end].nr);
1052     }
1053 }
1054
1055 /* enum for chain separation */
1056 enum ChainSepType {
1057     enChainSep_id_or_ter, enChainSep_id_and_ter, enChainSep_ter,
1058     enChainSep_id, enChainSep_interactive
1059 };
1060 static const char *ChainSepEnum[] = {"id_or_ter", "id_and_ter", "ter", "id", "interactive"};
1061 static const char *ChainSepInfoString[] =
1062 {
1063     "Splitting chemical chains based on TER records or chain id changing.\n",
1064     "Splitting chemical chains based on TER records and chain id changing.\n",
1065     "Splitting chemical chains based on TER records only (ignoring chain id).\n",
1066     "Splitting chemical chains based on changing chain id only (ignoring TER records).\n",
1067     "Splitting chemical chains interactively.\n"
1068 };
1069
1070 static void
1071 modify_chain_numbers(t_atoms *       pdba,
1072                      ChainSepType    enumChainSep)
1073 {
1074     int           i;
1075     char          old_prev_chainid;
1076     char          old_this_chainid;
1077     int           old_prev_chainnum;
1078     int           old_this_chainnum;
1079     t_resinfo    *ri;
1080     char          select[STRLEN];
1081     int           new_chainnum;
1082     int           this_atomnum;
1083     int           prev_atomnum;
1084     const char *  prev_atomname;
1085     const char *  this_atomname;
1086     const char *  prev_resname;
1087     const char *  this_resname;
1088     int           prev_resnum;
1089     int           this_resnum;
1090     char          prev_chainid;
1091     char          this_chainid;
1092
1093     /* The default chain enumeration is based on TER records only */
1094     printf("%s", ChainSepInfoString[enumChainSep]);
1095
1096     old_prev_chainid  = '?';
1097     old_prev_chainnum = -1;
1098     new_chainnum      = -1;
1099
1100     this_atomname       = nullptr;
1101     this_atomnum        = -1;
1102     this_resname        = nullptr;
1103     this_resnum         = -1;
1104     this_chainid        = '?';
1105
1106     for (i = 0; i < pdba->nres; i++)
1107     {
1108         ri                 = &pdba->resinfo[i];
1109         old_this_chainid   = ri->chainid;
1110         old_this_chainnum  = ri->chainnum;
1111
1112         prev_atomname      = this_atomname;
1113         prev_atomnum       = this_atomnum;
1114         prev_resname       = this_resname;
1115         prev_resnum        = this_resnum;
1116         prev_chainid       = this_chainid;
1117
1118         this_atomname      = *(pdba->atomname[i]);
1119         this_atomnum       = (pdba->pdbinfo != nullptr) ? pdba->pdbinfo[i].atomnr : i+1;
1120         this_resname       = *ri->name;
1121         this_resnum        = ri->nr;
1122         this_chainid       = ri->chainid;
1123
1124         switch (enumChainSep)
1125         {
1126             case enChainSep_id_or_ter:
1127                 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1128                 {
1129                     new_chainnum++;
1130                 }
1131                 break;
1132
1133             case enChainSep_id_and_ter:
1134                 if (old_this_chainid != old_prev_chainid && old_this_chainnum != old_prev_chainnum)
1135                 {
1136                     new_chainnum++;
1137                 }
1138                 break;
1139
1140             case enChainSep_id:
1141                 if (old_this_chainid != old_prev_chainid)
1142                 {
1143                     new_chainnum++;
1144                 }
1145                 break;
1146
1147             case enChainSep_ter:
1148                 if (old_this_chainnum != old_prev_chainnum)
1149                 {
1150                     new_chainnum++;
1151                 }
1152                 break;
1153             case enChainSep_interactive:
1154                 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1155                 {
1156                     if (i > 0)
1157                     {
1158                         printf("Split the chain (and introduce termini) between residue %s%d (chain id '%c', atom %d %s)\
1159 \n"
1160                                "and residue %s%d (chain id '%c', atom %d %s) ? [n/y]\n",
1161                                prev_resname, prev_resnum, prev_chainid, prev_atomnum, prev_atomname,
1162                                this_resname, this_resnum, this_chainid, this_atomnum, this_atomname);
1163
1164                         if (nullptr == fgets(select, STRLEN-1, stdin))
1165                         {
1166                             gmx_fatal(FARGS, "Error reading from stdin");
1167                         }
1168                     }
1169                     if (i == 0 || select[0] == 'y')
1170                     {
1171                         new_chainnum++;
1172                     }
1173                 }
1174                 break;
1175             default:
1176                 gmx_fatal(FARGS, "Internal inconsistency - this shouldn't happen...");
1177         }
1178         old_prev_chainid  = old_this_chainid;
1179         old_prev_chainnum = old_this_chainnum;
1180
1181         ri->chainnum = new_chainnum;
1182     }
1183 }
1184
1185 typedef struct {
1186     char  chainid;
1187     char  chainnum;
1188     int   start;
1189     int   natom;
1190     bool  bAllWat;
1191     int   nterpairs;
1192     int  *chainstart;
1193 } t_pdbchain;
1194
1195 typedef struct {
1196     char          chainid;
1197     int           chainnum;
1198     bool          bAllWat;
1199     int           nterpairs;
1200     int          *chainstart;
1201     t_hackblock **ntdb;
1202     t_hackblock **ctdb;
1203     int          *r_start;
1204     int          *r_end;
1205     t_atoms      *pdba;
1206     rvec         *x;
1207 } t_chain;
1208
1209 // TODO make all enums into scoped enums
1210 /* enum for vsites */
1211 enum VSitesType {
1212     enVSites_none, enVSites_hydrogens, enVSites_aromatics
1213 };
1214 static const char *VSitesEnum[] = {"none", "hydrogens", "aromatics"};
1215
1216 /* enum for water model */
1217 enum WaterType {
1218     enWater_select, enWater_none, enWater_spc, enWater_spce,
1219     enWater_tip3p, enWater_tip4p, enWater_tip5p, enWater_tips3p
1220 };
1221 static const char *WaterEnum[] = {
1222     "select", "none", "spc", "spce",
1223     "tip3p", "tip4p", "tip5p", "tips3p"
1224 };
1225
1226 /* enum for merge */
1227 enum MergeType {
1228     enMerge_no, enMerge_all, enMerge_interactive
1229 };
1230 static const char *MergeEnum[] = {"no", "all", "interactive"};
1231
1232 namespace gmx
1233 {
1234
1235 namespace
1236 {
1237
1238 class pdb2gmx : public ICommandLineOptionsModule
1239 {
1240     public:
1241         pdb2gmx() :
1242             bVsites_(FALSE), bPrevWat_(FALSE), bVsiteAromatics_(FALSE),
1243             enumChainSep_(enChainSep_id_or_ter),
1244             enumVSites_(enVSites_none),
1245             enumWater_(enWater_select),
1246             enumMerge_(enMerge_no),
1247             itp_file_(nullptr),
1248             nincl_(0),
1249             nmol_(0),
1250             incls_(nullptr),
1251             mols_(nullptr),
1252             mHmult_(0)
1253         {
1254         }
1255
1256         // From ICommandLineOptionsModule
1257         void init(CommandLineModuleSettings * /*settings*/) override
1258         {
1259         }
1260
1261         void initOptions(IOptionsContainer                 *options,
1262                          ICommandLineOptionsModuleSettings *settings) override;
1263
1264         void optionsFinished() override;
1265
1266         int run() override;
1267
1268     private:
1269         bool               bNewRTP_;
1270         bool               bInter_;
1271         bool               bCysMan_;
1272         bool               bLysMan_;
1273         bool               bAspMan_;
1274         bool               bGluMan_;
1275         bool               bHisMan_;
1276         bool               bGlnMan_;
1277         bool               bArgMan_;
1278         bool               bTerMan_;
1279         bool               bUnA_;
1280         bool               bHeavyH_;
1281         bool               bSort_;
1282         bool               bAllowMissing_;
1283         bool               bRemoveH_;
1284         bool               bDeuterate_;
1285         bool               bVerbose_;
1286         bool               bChargeGroups_;
1287         bool               bCmap_;
1288         bool               bRenumRes_;
1289         bool               bRTPresname_;
1290         bool               bIndexSet_;
1291         bool               bOutputSet_;
1292         bool               bVsites_;
1293         bool               bWat_;
1294         bool               bPrevWat_;
1295         bool               bITP_;
1296         bool               bVsiteAromatics_;
1297         real               angle_;
1298         real               distance_;
1299         real               posre_fc_;
1300         real               long_bond_dist_;
1301         real               short_bond_dist_;
1302
1303         std::string        indexOutputFile_;
1304         std::string        outputFile_;
1305         std::string        topologyFile_;
1306         std::string        includeTopologyFile_;
1307         std::string        outputConfFile_;
1308         std::string        inputConfFile_;
1309         std::string        outFile_;
1310         std::string        ff_;
1311
1312         ChainSepType       enumChainSep_;
1313         VSitesType         enumVSites_;
1314         WaterType          enumWater_;
1315         MergeType          enumMerge_;
1316
1317         FILE              *itp_file_;
1318         char               forcefield_[STRLEN];
1319         char               ffdir_[STRLEN];
1320         char              *ffname_;
1321         char              *watermodel_;
1322         int                nincl_;
1323         int                nmol_;
1324         char             **incls_;
1325         t_mols            *mols_;
1326         real               mHmult_;
1327 };
1328
1329 void pdb2gmx::initOptions(IOptionsContainer                 *options,
1330                           ICommandLineOptionsModuleSettings *settings)
1331 {
1332     const char *desc[] = {
1333         "[THISMODULE] reads a [REF].pdb[ref] (or [REF].gro[ref]) file, reads",
1334         "some database files, adds hydrogens to the molecules and generates",
1335         "coordinates in GROMACS (GROMOS), or optionally [REF].pdb[ref], format",
1336         "and a topology in GROMACS format.",
1337         "These files can subsequently be processed to generate a run input file.",
1338         "[PAR]",
1339         "[THISMODULE] will search for force fields by looking for",
1340         "a [TT]forcefield.itp[tt] file in subdirectories [TT]<forcefield>.ff[tt]",
1341         "of the current working directory and of the GROMACS library directory",
1342         "as inferred from the path of the binary or the [TT]GMXLIB[tt] environment",
1343         "variable.",
1344         "By default the forcefield selection is interactive,",
1345         "but you can use the [TT]-ff[tt] option to specify one of the short names",
1346         "in the list on the command line instead. In that case [THISMODULE] just looks",
1347         "for the corresponding [TT]<forcefield>.ff[tt] directory.",
1348         "[PAR]",
1349         "After choosing a force field, all files will be read only from",
1350         "the corresponding force field directory.",
1351         "If you want to modify or add a residue types, you can copy the force",
1352         "field directory from the GROMACS library directory to your current",
1353         "working directory. If you want to add new protein residue types,",
1354         "you will need to modify [TT]residuetypes.dat[tt] in the library directory",
1355         "or copy the whole library directory to a local directory and set",
1356         "the environment variable [TT]GMXLIB[tt] to the name of that directory.",
1357         "Check Chapter 5 of the manual for more information about file formats.",
1358         "[PAR]",
1359
1360         "Note that a [REF].pdb[ref] file is nothing more than a file format, and it",
1361         "need not necessarily contain a protein structure. Every kind of",
1362         "molecule for which there is support in the database can be converted.",
1363         "If there is no support in the database, you can add it yourself.[PAR]",
1364
1365         "The program has limited intelligence, it reads a number of database",
1366         "files, that allow it to make special bonds (Cys-Cys, Heme-His, etc.),",
1367         "if necessary this can be done manually. The program can prompt the",
1368         "user to select which kind of LYS, ASP, GLU, CYS or HIS residue is",
1369         "desired. For Lys the choice is between neutral (two protons on NZ) or",
1370         "protonated (three protons, default), for Asp and Glu unprotonated",
1371         "(default) or protonated, for His the proton can be either on ND1,",
1372         "on NE2 or on both. By default these selections are done automatically.",
1373         "For His, this is based on an optimal hydrogen bonding",
1374         "conformation. Hydrogen bonds are defined based on a simple geometric",
1375         "criterion, specified by the maximum hydrogen-donor-acceptor angle",
1376         "and donor-acceptor distance, which are set by [TT]-angle[tt] and",
1377         "[TT]-dist[tt] respectively.[PAR]",
1378
1379         "The protonation state of N- and C-termini can be chosen interactively",
1380         "with the [TT]-ter[tt] flag.  Default termini are ionized (NH3+ and COO-),",
1381         "respectively.  Some force fields support zwitterionic forms for chains of",
1382         "one residue, but for polypeptides these options should NOT be selected.",
1383         "The AMBER force fields have unique forms for the terminal residues,",
1384         "and these are incompatible with the [TT]-ter[tt] mechanism. You need",
1385         "to prefix your N- or C-terminal residue names with \"N\" or \"C\"",
1386         "respectively to use these forms, making sure you preserve the format",
1387         "of the coordinate file. Alternatively, use named terminating residues",
1388         "(e.g. ACE, NME).[PAR]",
1389
1390         "The separation of chains is not entirely trivial since the markup",
1391         "in user-generated PDB files frequently varies and sometimes it",
1392         "is desirable to merge entries across a TER record, for instance",
1393         "if you want a disulfide bridge or distance restraints between",
1394         "two protein chains or if you have a HEME group bound to a protein.",
1395         "In such cases multiple chains should be contained in a single",
1396         "[TT]moleculetype[tt] definition.",
1397         "To handle this, [THISMODULE] uses two separate options.",
1398         "First, [TT]-chainsep[tt] allows you to choose when a new chemical chain should",
1399         "start, and termini added when applicable. This can be done based on the",
1400         "existence of TER records, when the chain id changes, or combinations of either",
1401         "or both of these. You can also do the selection fully interactively.",
1402         "In addition, there is a [TT]-merge[tt] option that controls how multiple chains",
1403         "are merged into one moleculetype, after adding all the chemical termini (or not).",
1404         "This can be turned off (no merging), all non-water chains can be merged into a",
1405         "single molecule, or the selection can be done interactively.[PAR]",
1406
1407         "[THISMODULE] will also check the occupancy field of the [REF].pdb[ref] file.",
1408         "If any of the occupancies are not one, indicating that the atom is",
1409         "not resolved well in the structure, a warning message is issued.",
1410         "When a [REF].pdb[ref] file does not originate from an X-ray structure determination",
1411         "all occupancy fields may be zero. Either way, it is up to the user",
1412         "to verify the correctness of the input data (read the article!).[PAR]",
1413
1414         "During processing the atoms will be reordered according to GROMACS",
1415         "conventions. With [TT]-n[tt] an index file can be generated that",
1416         "contains one group reordered in the same way. This allows you to",
1417         "convert a GROMOS trajectory and coordinate file to GROMOS. There is",
1418         "one limitation: reordering is done after the hydrogens are stripped",
1419         "from the input and before new hydrogens are added. This means that",
1420         "you should not use [TT]-ignh[tt].[PAR]",
1421
1422         "The [REF].gro[ref] and [TT].g96[tt] file formats do not support chain",
1423         "identifiers. Therefore it is useful to enter a [REF].pdb[ref] file name at",
1424         "the [TT]-o[tt] option when you want to convert a multi-chain [REF].pdb[ref] file.",
1425         "[PAR]",
1426
1427         "The option [TT]-vsite[tt] removes hydrogen and fast improper dihedral",
1428         "motions. Angular and out-of-plane motions can be removed by changing",
1429         "hydrogens into virtual sites and fixing angles, which fixes their",
1430         "position relative to neighboring atoms. Additionally, all atoms in the",
1431         "aromatic rings of the standard amino acids (i.e. PHE, TRP, TYR and HIS)",
1432         "can be converted into virtual sites, eliminating the fast improper dihedral",
1433         "fluctuations in these rings (but this feature is deprecated).",
1434         "[BB]Note[bb] that in this case all other hydrogen",
1435         "atoms are also converted to virtual sites. The mass of all atoms that are",
1436         "converted into virtual sites, is added to the heavy atoms.[PAR]",
1437         "Also slowing down of dihedral motion can be done with [TT]-heavyh[tt]",
1438         "done by increasing the hydrogen-mass by a factor of 4. This is also",
1439         "done for water hydrogens to slow down the rotational motion of water.",
1440         "The increase in mass of the hydrogens is subtracted from the bonded",
1441         "(heavy) atom so that the total mass of the system remains the same."
1442     };
1443
1444     settings->setHelpText(desc);
1445
1446     options->addOption(BooleanOption("newrtp")
1447                            .store(&bNewRTP_).defaultValue(false).hidden()
1448                            .description("Write the residue database in new format to [TT]new.rtp[tt]"));
1449     options->addOption(RealOption("lb")
1450                            .store(&long_bond_dist_).defaultValue(0.25).hidden()
1451                            .description("Long bond warning distance"));
1452     options->addOption(RealOption("sb")
1453                            .store(&short_bond_dist_).defaultValue(0.05).hidden()
1454                            .description("Short bond warning distance"));
1455     options->addOption(EnumOption<ChainSepType>("chainsep").enumValue(ChainSepEnum)
1456                            .store(&enumChainSep_)
1457                            .description("Condition in PDB files when a new chain should be started (adding termini)"));
1458     options->addOption(EnumOption<MergeType>("merge").enumValue(MergeEnum)
1459                            .store(&enumMerge_)
1460                            .description("Merge multiple chains into a single [moleculetype]"));
1461     options->addOption(StringOption("ff")
1462                            .store(&ff_).defaultValue("select")
1463                            .description("Force field, interactive by default. Use [TT]-h[tt] for information."));
1464     options->addOption(EnumOption<WaterType>("water")
1465                            .store(&enumWater_).enumValue(WaterEnum)
1466                            .description("Water model to use"));
1467     options->addOption(BooleanOption("inter")
1468                            .store(&bInter_).defaultValue(false)
1469                            .description("Set the next 8 options to interactive"));
1470     options->addOption(BooleanOption("ss")
1471                            .store(&bCysMan_).defaultValue(false)
1472                            .description("Interactive SS bridge selection"));
1473     options->addOption(BooleanOption("ter")
1474                            .store(&bTerMan_).defaultValue(false)
1475                            .description("Interactive termini selection, instead of charged (default)"));
1476     options->addOption(BooleanOption("lys")
1477                            .store(&bLysMan_).defaultValue(false)
1478                            .description("Interactive lysine selection, instead of charged"));
1479     options->addOption(BooleanOption("arg")
1480                            .store(&bArgMan_).defaultValue(false)
1481                            .description("Interactive arginine selection, instead of charged"));
1482     options->addOption(BooleanOption("asp")
1483                            .store(&bAspMan_).defaultValue(false)
1484                            .description("Interactive aspartic acid selection, instead of charged"));
1485     options->addOption(BooleanOption("glu")
1486                            .store(&bGluMan_).defaultValue(false)
1487                            .description("Interactive glutamic acid selection, instead of charged"));
1488     options->addOption(BooleanOption("gln")
1489                            .store(&bGlnMan_).defaultValue(false)
1490                            .description("Interactive glutamine selection, instead of charged"));
1491     options->addOption(BooleanOption("his")
1492                            .store(&bHisMan_).defaultValue(false)
1493                            .description("Interactive histidine selection, instead of checking H-bonds"));
1494     options->addOption(RealOption("angle")
1495                            .store(&angle_).defaultValue(135.0)
1496                            .description("Minimum hydrogen-donor-acceptor angle for a H-bond (degrees)"));
1497     options->addOption(RealOption("dist")
1498                            .store(&distance_).defaultValue(0.3)
1499                            .description("Maximum donor-acceptor distance for a H-bond (nm)"));
1500     options->addOption(BooleanOption("una")
1501                            .store(&bUnA_).defaultValue(false)
1502                            .description("Select aromatic rings with united CH atoms on phenylalanine, tryptophane and tyrosine"));
1503     options->addOption(BooleanOption("sort")
1504                            .store(&bSort_).defaultValue(true).hidden()
1505                            .description("Sort the residues according to database, turning this off is dangerous as charge groups might be broken in parts"));
1506     options->addOption(BooleanOption("ignh")
1507                            .store(&bRemoveH_).defaultValue(false)
1508                            .description("Ignore hydrogen atoms that are in the coordinate file"));
1509     options->addOption(BooleanOption("missing")
1510                            .store(&bAllowMissing_).defaultValue(false)
1511                            .description("Continue when atoms are missing and bonds cannot be made, dangerous"));
1512     options->addOption(BooleanOption("v")
1513                            .store(&bVerbose_).defaultValue(false)
1514                            .description("Be slightly more verbose in messages"));
1515     options->addOption(RealOption("posrefc")
1516                            .store(&posre_fc_).defaultValue(1000)
1517                            .description("Force constant for position restraints"));
1518     options->addOption(EnumOption<VSitesType>("vsite")
1519                            .store(&enumVSites_).enumValue(VSitesEnum)
1520                            .description("Convert atoms to virtual sites"));
1521     options->addOption(BooleanOption("heavyh")
1522                            .store(&bHeavyH_).defaultValue(false)
1523                            .description("Make hydrogen atoms heavy"));
1524     options->addOption(BooleanOption("deuterate")
1525                            .store(&bDeuterate_).defaultValue(false)
1526                            .description("Change the mass of hydrogens to 2 amu"));
1527     options->addOption(BooleanOption("chargegrp")
1528                            .store(&bChargeGroups_).defaultValue(true)
1529                            .description("Use charge groups in the [REF].rtp[ref] file"));
1530     options->addOption(BooleanOption("cmap")
1531                            .store(&bCmap_).defaultValue(true)
1532                            .description("Use cmap torsions (if enabled in the [REF].rtp[ref] file)"));
1533     options->addOption(BooleanOption("renum")
1534                            .store(&bRenumRes_).defaultValue(false)
1535                            .description("Renumber the residues consecutively in the output"));
1536     options->addOption(BooleanOption("rtpres")
1537                            .store(&bRTPresname_).defaultValue(false)
1538                            .description("Use [REF].rtp[ref] entry names as residue names"));
1539     options->addOption(FileNameOption("f")
1540                            .legacyType(efSTX).inputFile()
1541                            .store(&inputConfFile_).required()
1542                            .defaultBasename("protein").defaultType(efPDB)
1543                            .description("Structure file"));
1544     options->addOption(FileNameOption("o")
1545                            .legacyType(efSTO).outputFile()
1546                            .store(&outputConfFile_).required()
1547                            .defaultBasename("conf")
1548                            .description("Structure file"));
1549     options->addOption(FileNameOption("p")
1550                            .legacyType(efTOP).outputFile()
1551                            .store(&topologyFile_).required()
1552                            .defaultBasename("topol")
1553                            .description("Topology file"));
1554     options->addOption(FileNameOption("i")
1555                            .legacyType(efITP).outputFile()
1556                            .store(&includeTopologyFile_).required()
1557                            .defaultBasename("posre")
1558                            .description("Include file for topology"));
1559     options->addOption(FileNameOption("n")
1560                            .legacyType(efNDX).outputFile()
1561                            .store(&indexOutputFile_).storeIsSet(&bIndexSet_)
1562                            .defaultBasename("index")
1563                            .description("Index file"));
1564     options->addOption(FileNameOption("q")
1565                            .legacyType(efSTO).outputFile()
1566                            .store(&outFile_).storeIsSet(&bOutputSet_)
1567                            .defaultBasename("clean").defaultType(efPDB)
1568                            .description("Structure file"));
1569 }
1570
1571 void pdb2gmx::optionsFinished()
1572 {
1573     if (inputConfFile_.empty())
1574     {
1575         GMX_THROW(InconsistentInputError("You must supply an input file"));
1576     }
1577     if (bInter_)
1578     {
1579         /* if anything changes here, also change description of -inter */
1580         bCysMan_ = true;
1581         bTerMan_ = true;
1582         bLysMan_ = true;
1583         bArgMan_ = true;
1584         bAspMan_ = true;
1585         bGluMan_ = true;
1586         bGlnMan_ = true;
1587         bHisMan_ = true;
1588     }
1589
1590     if (bHeavyH_)
1591     {
1592         mHmult_ = 4.0;
1593     }
1594     else if (bDeuterate_)
1595     {
1596         mHmult_ = 2.0;
1597     }
1598     else
1599     {
1600         mHmult_ = 1.0;
1601     }
1602
1603     /* Force field selection, interactive or direct */
1604     choose_ff(strcmp(ff_.c_str(), "select") == 0 ? nullptr : ff_.c_str(),
1605               forcefield_, sizeof(forcefield_),
1606               ffdir_, sizeof(ffdir_));
1607
1608     if (strlen(forcefield_) > 0)
1609     {
1610         ffname_    = forcefield_;
1611         ffname_[0] = std::toupper(ffname_[0]);
1612     }
1613     else
1614     {
1615         gmx_fatal(FARGS, "Empty forcefield string");
1616     }
1617 }
1618
1619 int pdb2gmx::run()
1620 {
1621     char         select[STRLEN];
1622     int          nssbonds;
1623     t_ssbond    *ssbonds;
1624     t_hackblock *hb_chain;
1625     t_restp     *restp_chain;
1626
1627     int          this_atomnum;
1628     int          prev_atomnum;
1629     const char  *prev_atomname;
1630     const char  *this_atomname;
1631     const char  *prev_resname;
1632     const char  *this_resname;
1633     int          prev_resnum;
1634     int          this_resnum;
1635     char         prev_chainid;
1636     char         this_chainid;
1637     int          prev_chainnumber;
1638     int          this_chainnumber;
1639     int          this_chainstart;
1640     int          prev_chainstart;
1641
1642     printf("\nUsing the %s force field in directory %s\n\n",
1643            ffname_, ffdir_);
1644
1645     choose_watermodel(WaterEnum[enumWater_], ffdir_, &watermodel_);
1646
1647     switch (enumVSites_)
1648     {
1649         case enVSites_none:
1650             bVsites_         = false;
1651             bVsiteAromatics_ = false;
1652             break;
1653         case enVSites_hydrogens:
1654             bVsites_         = true;
1655             bVsiteAromatics_ = false;
1656             break;
1657         case enVSites_aromatics:
1658             bVsites_         = true;
1659             bVsiteAromatics_ = true;
1660             break;
1661         default:
1662             gmx_fatal(FARGS, "Internal inconsistency: VSitesEnum='%s'", VSitesEnum[enumVSites_]);
1663     }       /* end switch */
1664
1665     /* Open the symbol table */
1666     t_symtab symtab;
1667     open_symtab(&symtab);
1668
1669     /* Residue type database */
1670     ResidueType rt;
1671
1672     /* Read residue renaming database(s), if present */
1673     std::vector<std::string> rrn = fflib_search_file_end(ffdir_, ".r2b", FALSE);
1674
1675     int                      nrtprename = 0;
1676     rtprename_t             *rtprename  = nullptr;
1677     for (const auto &filename : rrn)
1678     {
1679         printf("going to rename %s\n", filename.c_str());
1680         FILE *fp = fflib_open(filename);
1681         read_rtprename(filename.c_str(), fp, &nrtprename, &rtprename);
1682         gmx_ffclose(fp);
1683     }
1684
1685     /* Add all alternative names from the residue renaming database to the list
1686        of recognized amino/nucleic acids. */
1687     for (int i = 0; i < nrtprename; i++)
1688     {
1689         /* Only add names if the 'standard' gromacs/iupac base name was found */
1690
1691         /* TODO this should be changed with gmx::optional so that we only need
1692          * to search rt once.
1693          */
1694         if (rt.nameIndexedInResidueTypes(rtprename[i].gmx))
1695         {
1696             std::string restype = rt.typeNameForIndexedResidue(rtprename[i].gmx);
1697             rt.addResidue(rtprename[i].main, restype);
1698             rt.addResidue(rtprename[i].nter, restype);
1699             rt.addResidue(rtprename[i].cter, restype);
1700             rt.addResidue(rtprename[i].bter, restype);
1701         }
1702     }
1703
1704     matrix      box;
1705     const char *watres;
1706     clear_mat(box);
1707     if (watermodel_ != nullptr && (strstr(watermodel_, "4p") ||
1708                                    strstr(watermodel_, "4P")))
1709     {
1710         watres = "HO4";
1711     }
1712     else if (watermodel_ != nullptr && (strstr(watermodel_, "5p") ||
1713                                         strstr(watermodel_, "5P")))
1714     {
1715         watres = "HO5";
1716     }
1717     else
1718     {
1719         watres = "HOH";
1720     }
1721
1722     AtomProperties aps;
1723     char          *title;
1724     int            ePBC;
1725     t_atoms        pdba_all;
1726     rvec          *pdbx;
1727     int            natom = read_pdball(inputConfFile_.c_str(), bOutputSet_, outFile_.c_str(),
1728                                        &title, &pdba_all, &pdbx, &ePBC, box, bRemoveH_,
1729                                        &symtab, &rt, watres, &aps, bVerbose_);
1730
1731     if (natom == 0)
1732     {
1733         std::string message = formatString("No atoms found in pdb file %s\n", inputConfFile_.c_str());
1734         GMX_THROW(InconsistentInputError(message));
1735     }
1736
1737     printf("Analyzing pdb file\n");
1738     int nwaterchain = 0;
1739
1740     modify_chain_numbers(&pdba_all, enumChainSep_);
1741
1742     int nchainmerges        = 0;
1743
1744     this_atomname       = nullptr;
1745     this_atomnum        = -1;
1746     this_resname        = nullptr;
1747     this_resnum         = -1;
1748     this_chainid        = '?';
1749     this_chainnumber    = -1;
1750     this_chainstart     = 0;
1751     /* Keep the compiler happy */
1752     prev_chainstart     = 0;
1753
1754     int         numChains = 0;
1755     int         maxch     = 16;
1756     t_pdbchain *pdb_ch;
1757     snew(pdb_ch, maxch);
1758
1759     t_resinfo *ri;
1760     bool       bMerged = false;
1761     for (int i = 0; (i < natom); i++)
1762     {
1763         ri = &pdba_all.resinfo[pdba_all.atom[i].resind];
1764
1765         /* TODO this should live in a helper object, and consolidate
1766            that with code in modify_chain_numbers */
1767         prev_atomname      = this_atomname;
1768         prev_atomnum       = this_atomnum;
1769         prev_resname       = this_resname;
1770         prev_resnum        = this_resnum;
1771         prev_chainid       = this_chainid;
1772         prev_chainnumber   = this_chainnumber;
1773         if (!bMerged)
1774         {
1775             prev_chainstart    = this_chainstart;
1776         }
1777
1778         this_atomname      = *pdba_all.atomname[i];
1779         this_atomnum       = (pdba_all.pdbinfo != nullptr) ? pdba_all.pdbinfo[i].atomnr : i+1;
1780         this_resname       = *ri->name;
1781         this_resnum        = ri->nr;
1782         this_chainid       = ri->chainid;
1783         this_chainnumber   = ri->chainnum;
1784
1785         bWat_ = gmx::equalCaseInsensitive(*ri->name, watres);
1786
1787         if ((i == 0) || (this_chainnumber != prev_chainnumber) || (bWat_ != bPrevWat_))
1788         {
1789             this_chainstart = pdba_all.atom[i].resind;
1790             bMerged         = false;
1791             if (i > 0 && !bWat_)
1792             {
1793                 if (!strncmp(MergeEnum[enumMerge_], "int", 3))
1794                 {
1795                     printf("Merge chain ending with residue %s%d (chain id '%c', atom %d %s) and chain starting with\n"
1796                            "residue %s%d (chain id '%c', atom %d %s) into a single moleculetype (keeping termini)? [n/y]\n",
1797                            prev_resname, prev_resnum, prev_chainid, prev_atomnum, prev_atomname,
1798                            this_resname, this_resnum, this_chainid, this_atomnum, this_atomname);
1799
1800                     if (nullptr == fgets(select, STRLEN-1, stdin))
1801                     {
1802                         gmx_fatal(FARGS, "Error reading from stdin");
1803                     }
1804                     bMerged = (select[0] == 'y');
1805                 }
1806                 else if (!strncmp(MergeEnum[enumMerge_], "all", 3))
1807                 {
1808                     bMerged = true;
1809                 }
1810             }
1811
1812             if (bMerged)
1813             {
1814                 pdb_ch[numChains-1].chainstart[pdb_ch[numChains-1].nterpairs] =
1815                     pdba_all.atom[i].resind - prev_chainstart;
1816                 pdb_ch[numChains-1].nterpairs++;
1817                 srenew(pdb_ch[numChains-1].chainstart, pdb_ch[numChains-1].nterpairs+1);
1818                 nchainmerges++;
1819             }
1820             else
1821             {
1822                 /* set natom for previous chain */
1823                 if (numChains > 0)
1824                 {
1825                     pdb_ch[numChains-1].natom = i-pdb_ch[numChains-1].start;
1826                 }
1827                 if (bWat_)
1828                 {
1829                     nwaterchain++;
1830                     ri->chainid = ' ';
1831                 }
1832                 /* check if chain identifier was used before */
1833                 for (int j = 0; (j < numChains); j++)
1834                 {
1835                     if (pdb_ch[j].chainid != ' ' && pdb_ch[j].chainid == ri->chainid)
1836                     {
1837                         printf("WARNING: Chain identifier '%c' is used in two non-sequential blocks.\n"
1838                                "They will be treated as separate chains unless you reorder your file.\n",
1839                                ri->chainid);
1840                     }
1841                 }
1842                 // TODO This is too convoluted. Use a std::vector
1843                 if (numChains == maxch)
1844                 {
1845                     maxch += 16;
1846                     srenew(pdb_ch, maxch);
1847                 }
1848                 pdb_ch[numChains].chainid  = ri->chainid;
1849                 pdb_ch[numChains].chainnum = ri->chainnum;
1850                 pdb_ch[numChains].start    = i;
1851                 pdb_ch[numChains].bAllWat  = bWat_;
1852                 if (bWat_)
1853                 {
1854                     pdb_ch[numChains].nterpairs = 0;
1855                 }
1856                 else
1857                 {
1858                     pdb_ch[numChains].nterpairs = 1;
1859                 }
1860                 snew(pdb_ch[numChains].chainstart, pdb_ch[numChains].nterpairs+1);
1861                 /* modified [numChains] to [0] below */
1862                 pdb_ch[numChains].chainstart[0] = 0;
1863                 numChains++;
1864             }
1865         }
1866         bPrevWat_ = bWat_;
1867     }
1868     pdb_ch[numChains-1].natom = natom-pdb_ch[numChains-1].start;
1869
1870     /* set all the water blocks at the end of the chain */
1871     int *swap_index;
1872     snew(swap_index, numChains);
1873     int  j = 0;
1874     for (int i = 0; i < numChains; i++)
1875     {
1876         if (!pdb_ch[i].bAllWat)
1877         {
1878             swap_index[j] = i;
1879             j++;
1880         }
1881     }
1882     for (int i = 0; i < numChains; i++)
1883     {
1884         if (pdb_ch[i].bAllWat)
1885         {
1886             swap_index[j] = i;
1887             j++;
1888         }
1889     }
1890     if (nwaterchain > 1)
1891     {
1892         printf("Moved all the water blocks to the end\n");
1893     }
1894
1895     t_atoms *pdba;
1896     t_chain *chains;
1897     snew(chains, numChains);
1898     /* copy pdb data and x for all chains */
1899     for (int i = 0; (i < numChains); i++)
1900     {
1901         int si                   = swap_index[i];
1902         chains[i].chainid    = pdb_ch[si].chainid;
1903         chains[i].chainnum   = pdb_ch[si].chainnum;
1904         chains[i].bAllWat    = pdb_ch[si].bAllWat;
1905         chains[i].nterpairs  = pdb_ch[si].nterpairs;
1906         chains[i].chainstart = pdb_ch[si].chainstart;
1907         snew(chains[i].ntdb, pdb_ch[si].nterpairs);
1908         snew(chains[i].ctdb, pdb_ch[si].nterpairs);
1909         snew(chains[i].r_start, pdb_ch[si].nterpairs);
1910         snew(chains[i].r_end, pdb_ch[si].nterpairs);
1911
1912         snew(chains[i].pdba, 1);
1913         init_t_atoms(chains[i].pdba, pdb_ch[si].natom, true);
1914         snew(chains[i].x, chains[i].pdba->nr);
1915         for (j = 0; j < chains[i].pdba->nr; j++)
1916         {
1917             chains[i].pdba->atom[j] = pdba_all.atom[pdb_ch[si].start+j];
1918             snew(chains[i].pdba->atomname[j], 1);
1919             *chains[i].pdba->atomname[j] =
1920                 gmx_strdup(*pdba_all.atomname[pdb_ch[si].start+j]);
1921             chains[i].pdba->pdbinfo[j] = pdba_all.pdbinfo[pdb_ch[si].start+j];
1922             copy_rvec(pdbx[pdb_ch[si].start+j], chains[i].x[j]);
1923         }
1924         /* Re-index the residues assuming that the indices are continuous */
1925         int k                = chains[i].pdba->atom[0].resind;
1926         int nres             = chains[i].pdba->atom[chains[i].pdba->nr-1].resind - k + 1;
1927         chains[i].pdba->nres = nres;
1928         for (int j = 0; j < chains[i].pdba->nr; j++)
1929         {
1930             chains[i].pdba->atom[j].resind -= k;
1931         }
1932         srenew(chains[i].pdba->resinfo, nres);
1933         for (int j = 0; j < nres; j++)
1934         {
1935             chains[i].pdba->resinfo[j] = pdba_all.resinfo[k+j];
1936             snew(chains[i].pdba->resinfo[j].name, 1);
1937             *chains[i].pdba->resinfo[j].name = gmx_strdup(*pdba_all.resinfo[k+j].name);
1938             /* make all chain identifiers equal to that of the chain */
1939             chains[i].pdba->resinfo[j].chainid = pdb_ch[si].chainid;
1940         }
1941     }
1942
1943     if (nchainmerges > 0)
1944     {
1945         printf("\nMerged chains into joint molecule definitions at %d places.\n\n",
1946                nchainmerges);
1947     }
1948
1949     printf("There are %d chains and %d blocks of water and "
1950            "%d residues with %d atoms\n",
1951            numChains-nwaterchain, nwaterchain,
1952            pdba_all.nres, natom);
1953
1954     printf("\n  %5s  %4s %6s\n", "chain", "#res", "#atoms");
1955     for (int i = 0; (i < numChains); i++)
1956     {
1957         printf("  %d '%c' %5d %6d  %s\n",
1958                i+1, chains[i].chainid ? chains[i].chainid : '-',
1959                chains[i].pdba->nres, chains[i].pdba->nr,
1960                chains[i].bAllWat ? "(only water)" : "");
1961     }
1962     printf("\n");
1963
1964     check_occupancy(&pdba_all, inputConfFile_.c_str(), bVerbose_);
1965
1966     /* Read atomtypes... */
1967     gpp_atomtype *atype = read_atype(ffdir_, &symtab);
1968
1969     /* read residue database */
1970     printf("Reading residue database... (%s)\n", forcefield_);
1971     std::vector<std::string> rtpf  = fflib_search_file_end(ffdir_, ".rtp", true);
1972     int                      nrtp  = 0;
1973     t_restp                 *restp = nullptr;
1974     for (const auto &filename : rtpf)
1975     {
1976         read_resall(filename.c_str(), &nrtp, &restp, atype, &symtab, false);
1977     }
1978     if (bNewRTP_)
1979     {
1980         /* Not correct with multiple rtp input files with different bonded types */
1981         FILE *fp = gmx_fio_fopen("new.rtp", "w");
1982         print_resall(fp, nrtp, restp, atype);
1983         gmx_fio_fclose(fp);
1984     }
1985
1986     /* read hydrogen database */
1987     t_hackblock       *ah;
1988     int                nah = read_h_db(ffdir_, &ah);
1989
1990     /* Read Termini database... */
1991     int                ntdblist;
1992     t_hackblock       *ntdb;
1993     t_hackblock       *ctdb;
1994     t_hackblock      **tdblist;
1995     int                nNtdb = read_ter_db(ffdir_, 'n', &ntdb, atype);
1996     int                nCtdb = read_ter_db(ffdir_, 'c', &ctdb, atype);
1997
1998     FILE              *top_file = gmx_fio_fopen(topologyFile_.c_str(), "w");
1999
2000     print_top_header(top_file, topologyFile_.c_str(), FALSE, ffdir_, mHmult_);
2001
2002     t_chain *cc;
2003     rvec    *x;
2004     for (int chain = 0; (chain < numChains); chain++)
2005     {
2006         cc = &(chains[chain]);
2007
2008         /* set pdba, natom and nres to the current chain */
2009         pdba     = cc->pdba;
2010         x        = cc->x;
2011         natom    = cc->pdba->nr;
2012         int nres = cc->pdba->nres;
2013
2014         if (cc->chainid && ( cc->chainid != ' ' ) )
2015         {
2016             printf("Processing chain %d '%c' (%d atoms, %d residues)\n",
2017                    chain+1, cc->chainid, natom, nres);
2018         }
2019         else
2020         {
2021             printf("Processing chain %d (%d atoms, %d residues)\n",
2022                    chain+1, natom, nres);
2023         }
2024
2025         process_chain(pdba, x, bUnA_, bUnA_, bUnA_, bLysMan_, bAspMan_, bGluMan_,
2026                       bHisMan_, bArgMan_, bGlnMan_, angle_, distance_, &symtab,
2027                       nrtprename, rtprename);
2028
2029         cc->chainstart[cc->nterpairs] = pdba->nres;
2030         j = 0;
2031         for (int i = 0; i < cc->nterpairs; i++)
2032         {
2033             find_nc_ter(pdba, cc->chainstart[i], cc->chainstart[i+1],
2034                         &(cc->r_start[j]), &(cc->r_end[j]), &rt);
2035
2036             if (cc->r_start[j] >= 0 && cc->r_end[j] >= 0)
2037             {
2038                 j++;
2039             }
2040         }
2041         cc->nterpairs = j;
2042         if (cc->nterpairs == 0)
2043         {
2044             printf("Problem with chain definition, or missing terminal residues.\n"
2045                    "This chain does not appear to contain a recognized chain molecule.\n"
2046                    "If this is incorrect, you can edit residuetypes.dat to modify the behavior.\n");
2047         }
2048
2049         /* Check for disulfides and other special bonds */
2050         nssbonds = mk_specbonds(pdba, x, bCysMan_, &ssbonds, bVerbose_);
2051
2052         if (nrtprename > 0)
2053         {
2054             rename_resrtp(pdba, cc->nterpairs, cc->r_start, cc->r_end, nrtprename, rtprename,
2055                           &symtab, bVerbose_);
2056         }
2057
2058         for (int i = 0; i < cc->nterpairs; i++)
2059         {
2060
2061             /* Set termini.
2062              * We first apply a filter so we only have the
2063              * termini that can be applied to the residue in question
2064              * (or a generic terminus if no-residue specific is available).
2065              */
2066             /* First the N terminus */
2067             if (nNtdb > 0)
2068             {
2069                 tdblist = filter_ter(nNtdb, ntdb,
2070                                      *pdba->resinfo[cc->r_start[i]].name,
2071                                      &ntdblist);
2072                 if (ntdblist == 0)
2073                 {
2074                     printf("No suitable end (N or 5') terminus found in database - assuming this residue\n"
2075                            "is already in a terminus-specific form and skipping terminus selection.\n");
2076                     cc->ntdb[i] = nullptr;
2077                 }
2078                 else
2079                 {
2080                     if (bTerMan_ && ntdblist > 1)
2081                     {
2082                         sprintf(select, "Select start terminus type for %s-%d",
2083                                 *pdba->resinfo[cc->r_start[i]].name,
2084                                 pdba->resinfo[cc->r_start[i]].nr);
2085                         cc->ntdb[i] = choose_ter(ntdblist, tdblist, select);
2086                     }
2087                     else
2088                     {
2089                         cc->ntdb[i] = tdblist[0];
2090                     }
2091
2092                     printf("Start terminus %s-%d: %s\n",
2093                            *pdba->resinfo[cc->r_start[i]].name,
2094                            pdba->resinfo[cc->r_start[i]].nr,
2095                            (cc->ntdb[i])->name);
2096                     sfree(tdblist);
2097                 }
2098             }
2099             else
2100             {
2101                 cc->ntdb[i] = nullptr;
2102             }
2103
2104             /* And the C terminus */
2105             if (nCtdb > 0)
2106             {
2107                 tdblist = filter_ter(nCtdb, ctdb,
2108                                      *pdba->resinfo[cc->r_end[i]].name,
2109                                      &ntdblist);
2110                 if (ntdblist == 0)
2111                 {
2112                     printf("No suitable end (C or 3') terminus found in database - assuming this residue\n"
2113                            "is already in a terminus-specific form and skipping terminus selection.\n");
2114                     cc->ctdb[i] = nullptr;
2115                 }
2116                 else
2117                 {
2118                     if (bTerMan_ && ntdblist > 1)
2119                     {
2120                         sprintf(select, "Select end terminus type for %s-%d",
2121                                 *pdba->resinfo[cc->r_end[i]].name,
2122                                 pdba->resinfo[cc->r_end[i]].nr);
2123                         cc->ctdb[i] = choose_ter(ntdblist, tdblist, select);
2124                     }
2125                     else
2126                     {
2127                         cc->ctdb[i] = tdblist[0];
2128                     }
2129                     printf("End terminus %s-%d: %s\n",
2130                            *pdba->resinfo[cc->r_end[i]].name,
2131                            pdba->resinfo[cc->r_end[i]].nr,
2132                            (cc->ctdb[i])->name);
2133                     sfree(tdblist);
2134                 }
2135             }
2136             else
2137             {
2138                 cc->ctdb[i] = nullptr;
2139             }
2140         }
2141
2142         /* lookup hackblocks and rtp for all residues */
2143         get_hackblocks_rtp(&hb_chain, &restp_chain,
2144                            nrtp, restp, pdba->nres, pdba->resinfo,
2145                            cc->nterpairs, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end,
2146                            bAllowMissing_);
2147         /* ideally, now we would not need the rtp itself anymore, but do
2148            everything using the hb and restp arrays. Unfortunately, that
2149            requires some re-thinking of code in gen_vsite.c, which I won't
2150            do now :( AF 26-7-99 */
2151
2152         rename_atoms(nullptr, ffdir_,
2153                      pdba, &symtab, restp_chain, false, &rt, false, bVerbose_);
2154
2155         match_atomnames_with_rtp(restp_chain, hb_chain, pdba, x, bVerbose_);
2156
2157         if (bSort_)
2158         {
2159             t_blocka  *block;
2160             char     **gnames;
2161             block = new_blocka();
2162             snew(gnames, 1);
2163             sort_pdbatoms(restp_chain, natom, &pdba, &x, block, &gnames);
2164             remove_duplicate_atoms(pdba, x, bVerbose_);
2165             if (bIndexSet_)
2166             {
2167                 if (bRemoveH_)
2168                 {
2169                     fprintf(stderr, "WARNING: with the -remh option the generated "
2170                             "index file (%s) might be useless\n"
2171                             "(the index file is generated before hydrogens are added)",
2172                             indexOutputFile_.c_str());
2173                 }
2174                 write_index(indexOutputFile_.c_str(), block, gnames, false, 0);
2175             }
2176             for (int i = 0; i < block->nr; i++)
2177             {
2178                 sfree(gnames[i]);
2179             }
2180             sfree(gnames);
2181             done_blocka(block);
2182         }
2183         else
2184         {
2185             fprintf(stderr, "WARNING: "
2186                     "without sorting no check for duplicate atoms can be done\n");
2187         }
2188
2189         /* Generate Hydrogen atoms (and termini) in the sequence */
2190         printf("Generating any missing hydrogen atoms and/or adding termini.\n");
2191         add_h(&pdba, &x, nah, ah,
2192               cc->nterpairs, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end, bAllowMissing_,
2193               nullptr, nullptr, true, false);
2194         printf("Now there are %d residues with %d atoms\n",
2195                pdba->nres, pdba->nr);
2196
2197         /* make up molecule name(s) */
2198
2199         int         k = (cc->nterpairs > 0 && cc->r_start[0] >= 0) ? cc->r_start[0] : 0;
2200
2201         std::string restype = rt.typeNameForIndexedResidue(*pdba->resinfo[k].name);
2202
2203         std::string molname;
2204         std::string suffix;
2205         if (cc->bAllWat)
2206         {
2207             molname = "Water";
2208         }
2209         else
2210         {
2211             this_chainid = cc->chainid;
2212
2213             /* Add the chain id if we have one */
2214             if (this_chainid != ' ')
2215             {
2216                 suffix.append(formatString("_chain_%c", this_chainid));
2217             }
2218
2219             /* Check if there have been previous chains with the same id */
2220             int nid_used = 0;
2221             for (int k = 0; k < chain; k++)
2222             {
2223                 if (cc->chainid == chains[k].chainid)
2224                 {
2225                     nid_used++;
2226                 }
2227             }
2228             /* Add the number for this chain identifier if there are multiple copies */
2229             if (nid_used > 0)
2230             {
2231                 suffix.append(formatString("%d", nid_used+1));
2232             }
2233
2234             if (suffix.length() > 0)
2235             {
2236                 molname.append(restype);
2237                 molname.append(suffix);
2238             }
2239             else
2240             {
2241                 molname = restype;
2242             }
2243         }
2244         std::string itp_fn   = topologyFile_;;
2245         std::string posre_fn = includeTopologyFile_;
2246         if ((numChains-nwaterchain > 1) && !cc->bAllWat)
2247         {
2248             bITP_  = true;
2249             printf("Chain time...\n");
2250             //construct the itp file name
2251             itp_fn = stripSuffixIfPresent(itp_fn, ".top");
2252             itp_fn.append("_");
2253             itp_fn.append(molname);
2254             itp_fn.append(".itp");
2255             //now do the same for posre
2256             posre_fn = stripSuffixIfPresent(posre_fn, ".itp");
2257             posre_fn.append("_");
2258             posre_fn.append(molname);
2259             posre_fn.append(".itp");
2260             if (posre_fn == itp_fn)
2261             {
2262                 posre_fn = Path::concatenateBeforeExtension(posre_fn, "_pr");
2263             }
2264             nincl_++;
2265             srenew(incls_, nincl_);
2266             incls_[nincl_-1] = gmx_strdup(itp_fn.c_str());
2267             itp_file_        = gmx_fio_fopen(itp_fn.c_str(), "w");
2268         }
2269         else
2270         {
2271             bITP_ = false;
2272         }
2273
2274         srenew(mols_, nmol_+1);
2275         if (cc->bAllWat)
2276         {
2277             mols_[nmol_].name = gmx_strdup("SOL");
2278             mols_[nmol_].nr   = pdba->nres;
2279         }
2280         else
2281         {
2282             mols_[nmol_].name = gmx_strdup(molname.c_str());
2283             mols_[nmol_].nr   = 1;
2284         }
2285         nmol_++;
2286
2287         if (bITP_)
2288         {
2289             print_top_comment(itp_file_, itp_fn.c_str(), ffdir_, true);
2290         }
2291
2292         FILE *top_file2;
2293         if (cc->bAllWat)
2294         {
2295             top_file2 = nullptr;
2296         }
2297         else if (bITP_)
2298         {
2299             top_file2 = itp_file_;
2300         }
2301         else
2302         {
2303             top_file2 = top_file;
2304         }
2305
2306         pdb2top(top_file2, posre_fn.c_str(), molname.c_str(), pdba, &x, atype, &symtab,
2307                 nrtp, restp,
2308                 restp_chain, hb_chain,
2309                 bAllowMissing_,
2310                 bVsites_, bVsiteAromatics_, ffdir_,
2311                 mHmult_, nssbonds, ssbonds,
2312                 long_bond_dist_, short_bond_dist_, bDeuterate_, bChargeGroups_, bCmap_,
2313                 bRenumRes_, bRTPresname_);
2314
2315         if (!cc->bAllWat)
2316         {
2317             write_posres(posre_fn.c_str(), pdba, posre_fc_);
2318         }
2319
2320         if (bITP_)
2321         {
2322             gmx_fio_fclose(itp_file_);
2323         }
2324
2325         /* pdba and natom have been reassigned somewhere so: */
2326         cc->pdba = pdba;
2327         cc->x    = x;
2328
2329     }
2330
2331     if (watermodel_ == nullptr)
2332     {
2333         for (int chain = 0; chain < numChains; chain++)
2334         {
2335             if (chains[chain].bAllWat)
2336             {
2337                 auto message = formatString("You have chosen not to include a water model, "
2338                                             "but there is water in the input file. Select a "
2339                                             "water model or remove the water from your input file.");
2340                 GMX_THROW(InconsistentInputError(message));
2341             }
2342         }
2343     }
2344     else
2345     {
2346         std::string waterFile = formatString("%s%c%s.itp", ffdir_, DIR_SEPARATOR, watermodel_);
2347         if (!fflib_fexist(waterFile))
2348         {
2349             auto message = formatString("The topology file '%s' for the selected water "
2350                                         "model '%s' can not be found in the force field "
2351                                         "directory. Select a different water model.",
2352                                         waterFile.c_str(), watermodel_);
2353             GMX_THROW(InconsistentInputError(message));
2354         }
2355     }
2356
2357     print_top_mols(top_file, title, ffdir_, watermodel_, nincl_, incls_, nmol_, mols_);
2358     gmx_fio_fclose(top_file);
2359
2360     /* now merge all chains back together */
2361     natom     = 0;
2362     int nres  = 0;
2363     for (int i = 0; (i < numChains); i++)
2364     {
2365         natom += chains[i].pdba->nr;
2366         nres  += chains[i].pdba->nres;
2367     }
2368     t_atoms *atoms;
2369     snew(atoms, 1);
2370     init_t_atoms(atoms, natom, false);
2371     for (int i = 0; i < atoms->nres; i++)
2372     {
2373         sfree(atoms->resinfo[i].name);
2374     }
2375     sfree(atoms->resinfo);
2376     atoms->nres = nres;
2377     snew(atoms->resinfo, nres);
2378     snew(x, natom);
2379     int k = 0;
2380     int l = 0;
2381     for (int i = 0; (i < numChains); i++)
2382     {
2383         if (numChains > 1)
2384         {
2385             printf("Including chain %d in system: %d atoms %d residues\n",
2386                    i+1, chains[i].pdba->nr, chains[i].pdba->nres);
2387         }
2388         for (int j = 0; (j < chains[i].pdba->nr); j++)
2389         {
2390             atoms->atom[k]         = chains[i].pdba->atom[j];
2391             atoms->atom[k].resind += l; /* l is processed nr of residues */
2392             atoms->atomname[k]     = chains[i].pdba->atomname[j];
2393             atoms->resinfo[atoms->atom[k].resind].chainid = chains[i].chainid;
2394             copy_rvec(chains[i].x[j], x[k]);
2395             k++;
2396         }
2397         for (int j = 0; (j < chains[i].pdba->nres); j++)
2398         {
2399             atoms->resinfo[l] = chains[i].pdba->resinfo[j];
2400             if (bRTPresname_)
2401             {
2402                 atoms->resinfo[l].name = atoms->resinfo[l].rtp;
2403             }
2404             l++;
2405         }
2406     }
2407
2408     if (numChains > 1)
2409     {
2410         fprintf(stderr, "Now there are %d atoms and %d residues\n", k, l);
2411         print_sums(atoms, true);
2412     }
2413
2414     rvec box_space;
2415     fprintf(stderr, "\nWriting coordinate file...\n");
2416     clear_rvec(box_space);
2417     if (box[0][0] == 0)
2418     {
2419         make_new_box(atoms->nr, x, box, box_space, false);
2420     }
2421     write_sto_conf(outputConfFile_.c_str(), title, atoms, x, nullptr, ePBC, box);
2422
2423     printf("\t\t--------- PLEASE NOTE ------------\n");
2424     printf("You have successfully generated a topology from: %s.\n",
2425            inputConfFile_.c_str());
2426     if (watermodel_ != nullptr)
2427     {
2428         printf("The %s force field and the %s water model are used.\n",
2429                ffname_, watermodel_);
2430     }
2431     else
2432     {
2433         printf("The %s force field is used.\n",
2434                ffname_);
2435     }
2436     printf("\t\t--------- ETON ESAELP ------------\n");
2437
2438     return 0;
2439 }
2440
2441 }   // namespace
2442
2443 const char pdb2gmxInfo::name[]             = "pdb2gmx";
2444 const char pdb2gmxInfo::shortDescription[] =
2445     "Convert coordinate files to topology and FF-compliant coordinate files";
2446 ICommandLineOptionsModulePointer pdb2gmxInfo::create()
2447 {
2448     return std::make_unique<pdb2gmx>();
2449 }
2450
2451 } // namespace gmx