src/gromacs/gmxpreprocess/pdb2gmx.c

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