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 "pdb2gmx.h"
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <ctype.h>
  44 #include <stdlib.h>
  45 #include <string.h>
  46 #include <time.h>
  47
  48 #include "typedefs.h"
  49 #include "gromacs/fileio/gmxfio.h"
  50 #include "copyrite.h"
  51 #include "gromacs/utility/cstringutil.h"
  52 #include "gromacs/fileio/confio.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/utility/futil.h"
  55 #include "gromacs/fileio/pdbio.h"
  56 #include "toputil.h"
  57 #include "h_db.h"
  58 #include "pgutil.h"
  59 #include "resall.h"
  60 #include "pdb2top.h"
  61 #include "ter_db.h"
  62 #include "gromacs/gmxlib/conformation-utilities.h"
  63 #include "genhydro.h"
  64 #include "readinp.h"
  65 #include "index.h"
  66 #include "fflibutil.h"
  67 #include "macros.h"
  68
  69 #include "gromacs/commandline/pargs.h"
  70 #include "gromacs/fileio/strdb.h"
  71 #include "gromacs/topology/atomprop.h"
  72 #include "gromacs/topology/block.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, gmx_residuetype_t rt)
 839 {
 840     int         i;
 841     const char *p_startrestype;
 842     const char *p_restype;
 843     int         nstartwarn, nendwarn;
 844
 845     *r_start = -1;
 846     *r_end   = -1;
 847
 848     nstartwarn = 0;
 849     nendwarn   = 0;
 850
 851     /* Find the starting terminus (typially N or 5') */
 852     for (i = r0; i < r1 && *r_start == -1; i++)
 853     {
 854         gmx_residuetype_get_type(rt, *pdba->resinfo[i].name, &p_startrestype);
 855         if (!gmx_strcasecmp(p_startrestype, "Protein") || !gmx_strcasecmp(p_startrestype, "DNA") || !gmx_strcasecmp(p_startrestype, "RNA") )
 856         {
 857             printf("Identified residue %s%d as a starting terminus.\n", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
 858             *r_start = i;
 859         }
 860         else
 861         {
 862             if (nstartwarn < 5)
 863             {
 864                 printf("Warning: Starting residue %s%d in chain not identified as Protein/RNA/DNA.\n", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
 865             }
 866             if (nstartwarn == 5)
 867             {
 868                 printf("More than 5 unidentified residues at start of chain - disabling further warnings.\n");
 869             }
 870             nstartwarn++;
 871         }
 872     }
 873
 874     if (*r_start >= 0)
 875     {
 876         /* Go through the rest of the residues, check that they are the same class, and identify the ending terminus. */
 877         for (i = *r_start; i < r1; i++)
 878         {
 879             gmx_residuetype_get_type(rt, *pdba->resinfo[i].name, &p_restype);
 880             if (!gmx_strcasecmp(p_restype, p_startrestype) && nendwarn == 0)
 881             {
 882                 *r_end = i;
 883             }
 884             else
 885             {
 886                 if (nendwarn < 5)
 887                 {
 888                     printf("Warning: Residue %s%d in chain has different type (%s) from starting residue %s%d (%s).\n",
 889                            *pdba->resinfo[i].name, pdba->resinfo[i].nr, p_restype,
 890                            *pdba->resinfo[*r_start].name, pdba->resinfo[*r_start].nr, p_startrestype);
 891                 }
 892                 if (nendwarn == 5)
 893                 {
 894                     printf("More than 5 unidentified residues at end of chain - disabling further warnings.\n");
 895                 }
 896                 nendwarn++;
 897             }
 898         }
 899     }
 900
 901     if (*r_end >= 0)
 902     {
 903         printf("Identified residue %s%d as a ending terminus.\n", *pdba->resinfo[*r_end].name, pdba->resinfo[*r_end].nr);
 904     }
 905 }
 906
 907
 908 static void
 909 modify_chain_numbers(t_atoms *       pdba,
 910                      const char *    chainsep)
 911 {
 912     int           i;
 913     char          old_prev_chainid;
 914     char          old_this_chainid;
 915     int           old_prev_chainnum;
 916     int           old_this_chainnum;
 917     t_resinfo    *ri;
 918     char          select[STRLEN];
 919     int           new_chainnum;
 920     int           this_atomnum;
 921     int           prev_atomnum;
 922     const char *  prev_atomname;
 923     const char *  this_atomname;
 924     const char *  prev_resname;
 925     const char *  this_resname;
 926     int           prev_resnum;
 927     int           this_resnum;
 928     char          prev_chainid;
 929     char          this_chainid;
 930     int           prev_chainnumber;
 931     int           this_chainnumber;
 932
 933     enum
 934     {
 935         SPLIT_ID_OR_TER,
 936         SPLIT_ID_AND_TER,
 937         SPLIT_ID_ONLY,
 938         SPLIT_TER_ONLY,
 939         SPLIT_INTERACTIVE
 940     }
 941     splitting;
 942
 943     splitting = SPLIT_TER_ONLY; /* keep compiler happy */
 944
 945     /* Be a bit flexible to catch typos */
 946     if (!strncmp(chainsep, "id_o", 4))
 947     {
 948         /* For later interactive splitting we tentatively assign new chain numbers at either changing id or ter records */
 949         splitting = SPLIT_ID_OR_TER;
 950         printf("Splitting chemical chains based on TER records or chain id changing.\n");
 951     }
 952     else if (!strncmp(chainsep, "int", 3))
 953     {
 954         /* For later interactive splitting we tentatively assign new chain numbers at either changing id or ter records */
 955         splitting = SPLIT_INTERACTIVE;
 956         printf("Splitting chemical chains interactively.\n");
 957     }
 958     else if (!strncmp(chainsep, "id_a", 4))
 959     {
 960         splitting = SPLIT_ID_AND_TER;
 961         printf("Splitting chemical chains based on TER records and chain id changing.\n");
 962     }
 963     else if (strlen(chainsep) == 2 && !strncmp(chainsep, "id", 4))
 964     {
 965         splitting = SPLIT_ID_ONLY;
 966         printf("Splitting chemical chains based on changing chain id only (ignoring TER records).\n");
 967     }
 968     else if (chainsep[0] == 't')
 969     {
 970         splitting = SPLIT_TER_ONLY;
 971         printf("Splitting chemical chains based on TER records only (ignoring chain id).\n");
 972     }
 973     else
 974     {
 975         gmx_fatal(FARGS, "Unidentified setting for chain separation: %s\n", chainsep);
 976     }
 977
 978     /* The default chain enumeration is based on TER records only, which is reflected in chainnum below */
 979
 980     old_prev_chainid  = '?';
 981     old_prev_chainnum = -1;
 982     new_chainnum      = -1;
 983
 984     this_atomname       = NULL;
 985     this_atomnum        = -1;
 986     this_resname        = NULL;
 987     this_resnum         = -1;
 988     this_chainid        = '?';
 989     this_chainnumber    = -1;
 990
 991     for (i = 0; i < pdba->nres; i++)
 992     {
 993         ri                 = &pdba->resinfo[i];
 994         old_this_chainid   = ri->chainid;
 995         old_this_chainnum  = ri->chainnum;
 996
 997         prev_atomname      = this_atomname;
 998         prev_atomnum       = this_atomnum;
 999         prev_resname       = this_resname;
1000         prev_resnum        = this_resnum;
1001         prev_chainid       = this_chainid;
1002         prev_chainnumber   = this_chainnumber;
1003
1004         this_atomname      = *(pdba->atomname[i]);
1005         this_atomnum       = (pdba->pdbinfo != NULL) ? pdba->pdbinfo[i].atomnr : i+1;
1006         this_resname       = *ri->name;
1007         this_resnum        = ri->nr;
1008         this_chainid       = ri->chainid;
1009         this_chainnumber   = ri->chainnum;
1010
1011         switch (splitting)
1012         {
1013             case SPLIT_ID_OR_TER:
1014                 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1015                 {
1016                     new_chainnum++;
1017                 }
1018                 break;
1019
1020             case SPLIT_ID_AND_TER:
1021                 if (old_this_chainid != old_prev_chainid && old_this_chainnum != old_prev_chainnum)
1022                 {
1023                     new_chainnum++;
1024                 }
1025                 break;
1026
1027             case SPLIT_ID_ONLY:
1028                 if (old_this_chainid != old_prev_chainid)
1029                 {
1030                     new_chainnum++;
1031                 }
1032                 break;
1033
1034             case SPLIT_TER_ONLY:
1035                 if (old_this_chainnum != old_prev_chainnum)
1036                 {
1037                     new_chainnum++;
1038                 }
1039                 break;
1040             case SPLIT_INTERACTIVE:
1041                 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1042                 {
1043                     if (i > 0)
1044                     {
1045                         printf("Split the chain (and introduce termini) between residue %s%d (chain id '%c', atom %d %s)\n"
1046                                "and residue %s%d (chain id '%c', atom %d %s) ? [n/y]\n",
1047                                prev_resname, prev_resnum, prev_chainid, prev_atomnum, prev_atomname,
1048                                this_resname, this_resnum, this_chainid, this_atomnum, this_atomname);
1049
1050                         if (NULL == fgets(select, STRLEN-1, stdin))
1051                         {
1052                             gmx_fatal(FARGS, "Error reading from stdin");
1053                         }
1054                     }
1055                     if (i == 0 || select[0] == 'y')
1056                     {
1057                         new_chainnum++;
1058                     }
1059                 }
1060                 break;
1061             default:
1062                 gmx_fatal(FARGS, "Internal inconsistency - this shouldn't happen...");
1063                 break;
1064         }
1065         old_prev_chainid  = old_this_chainid;
1066         old_prev_chainnum = old_this_chainnum;
1067
1068         ri->chainnum = new_chainnum;
1069     }
1070 }
1071
1072
1073 typedef struct {
1074     char     chainid;
1075     char     chainnum;
1076     int      start;
1077     int      natom;
1078     gmx_bool bAllWat;
1079     int      nterpairs;
1080     int     *chainstart;
1081 } t_pdbchain;
1082
1083 typedef struct {
1084     char          chainid;
1085     int           chainnum;
1086     gmx_bool      bAllWat;
1087     int           nterpairs;
1088     int          *chainstart;
1089     t_hackblock **ntdb;
1090     t_hackblock **ctdb;
1091     int          *r_start;
1092     int          *r_end;
1093     t_atoms      *pdba;
1094     rvec         *x;
1095 } t_chain;
1096
1097 int gmx_pdb2gmx(int argc, char *argv[])
1098 {
1099     const char *desc[] = {
1100         "[THISMODULE] reads a [TT].pdb[tt] (or [TT].gro[tt]) file, reads",
1101         "some database files, adds hydrogens to the molecules and generates",
1102         "coordinates in GROMACS (GROMOS), or optionally [TT].pdb[tt], format",
1103         "and a topology in GROMACS format.",
1104         "These files can subsequently be processed to generate a run input file.",
1105         "[PAR]",
1106         "[THISMODULE] will search for force fields by looking for",
1107         "a [TT]forcefield.itp[tt] file in subdirectories [TT]<forcefield>.ff[tt]",
1108         "of the current working directory and of the GROMACS library directory",
1109         "as inferred from the path of the binary or the [TT]GMXLIB[tt] environment",
1110         "variable.",
1111         "By default the forcefield selection is interactive,",
1112         "but you can use the [TT]-ff[tt] option to specify one of the short names",
1113         "in the list on the command line instead. In that case [THISMODULE] just looks",
1114         "for the corresponding [TT]<forcefield>.ff[tt] directory.",
1115         "[PAR]",
1116         "After choosing a force field, all files will be read only from",
1117         "the corresponding force field directory.",
1118         "If you want to modify or add a residue types, you can copy the force",
1119         "field directory from the GROMACS library directory to your current",
1120         "working directory. If you want to add new protein residue types,",
1121         "you will need to modify [TT]residuetypes.dat[tt] in the library directory",
1122         "or copy the whole library directory to a local directory and set",
1123         "the environment variable [TT]GMXLIB[tt] to the name of that directory.",
1124         "Check Chapter 5 of the manual for more information about file formats.",
1125         "[PAR]",
1126
1127         "Note that a [TT].pdb[tt] file is nothing more than a file format, and it",
1128         "need not necessarily contain a protein structure. Every kind of",
1129         "molecule for which there is support in the database can be converted.",
1130         "If there is no support in the database, you can add it yourself.[PAR]",
1131
1132         "The program has limited intelligence, it reads a number of database",
1133         "files, that allow it to make special bonds (Cys-Cys, Heme-His, etc.),",
1134         "if necessary this can be done manually. The program can prompt the",
1135         "user to select which kind of LYS, ASP, GLU, CYS or HIS residue is",
1136         "desired. For Lys the choice is between neutral (two protons on NZ) or",
1137         "protonated (three protons, default), for Asp and Glu unprotonated",
1138         "(default) or protonated, for His the proton can be either on ND1,",
1139         "on NE2 or on both. By default these selections are done automatically.",
1140         "For His, this is based on an optimal hydrogen bonding",
1141         "conformation. Hydrogen bonds are defined based on a simple geometric",
1142         "criterion, specified by the maximum hydrogen-donor-acceptor angle",
1143         "and donor-acceptor distance, which are set by [TT]-angle[tt] and",
1144         "[TT]-dist[tt] respectively.[PAR]",
1145
1146         "The protonation state of N- and C-termini can be chosen interactively",
1147         "with the [TT]-ter[tt] flag.  Default termini are ionized (NH3+ and COO-),",
1148         "respectively.  Some force fields support zwitterionic forms for chains of",
1149         "one residue, but for polypeptides these options should NOT be selected.",
1150         "The AMBER force fields have unique forms for the terminal residues,",
1151         "and these are incompatible with the [TT]-ter[tt] mechanism. You need",
1152         "to prefix your N- or C-terminal residue names with \"N\" or \"C\"",
1153         "respectively to use these forms, making sure you preserve the format",
1154         "of the coordinate file. Alternatively, use named terminating residues",
1155         "(e.g. ACE, NME).[PAR]",
1156
1157         "The separation of chains is not entirely trivial since the markup",
1158         "in user-generated PDB files frequently varies and sometimes it",
1159         "is desirable to merge entries across a TER record, for instance",
1160         "if you want a disulfide bridge or distance restraints between",
1161         "two protein chains or if you have a HEME group bound to a protein.",
1162         "In such cases multiple chains should be contained in a single",
1163         "[TT]moleculetype[tt] definition.",
1164         "To handle this, [THISMODULE] uses two separate options.",
1165         "First, [TT]-chainsep[tt] allows you to choose when a new chemical chain should",
1166         "start, and termini added when applicable. This can be done based on the",
1167         "existence of TER records, when the chain id changes, or combinations of either",
1168         "or both of these. You can also do the selection fully interactively.",
1169         "In addition, there is a [TT]-merge[tt] option that controls how multiple chains",
1170         "are merged into one moleculetype, after adding all the chemical termini (or not).",
1171         "This can be turned off (no merging), all non-water chains can be merged into a",
1172         "single molecule, or the selection can be done interactively.[PAR]",
1173
1174         "[THISMODULE] will also check the occupancy field of the [TT].pdb[tt] file.",
1175         "If any of the occupancies are not one, indicating that the atom is",
1176         "not resolved well in the structure, a warning message is issued.",
1177         "When a [TT].pdb[tt] file does not originate from an X-ray structure determination",
1178         "all occupancy fields may be zero. Either way, it is up to the user",
1179         "to verify the correctness of the input data (read the article!).[PAR]",
1180
1181         "During processing the atoms will be reordered according to GROMACS",
1182         "conventions. With [TT]-n[tt] an index file can be generated that",
1183         "contains one group reordered in the same way. This allows you to",
1184         "convert a GROMOS trajectory and coordinate file to GROMOS. There is",
1185         "one limitation: reordering is done after the hydrogens are stripped",
1186         "from the input and before new hydrogens are added. This means that",
1187         "you should not use [TT]-ignh[tt].[PAR]",
1188
1189         "The [TT].gro[tt] and [TT].g96[tt] file formats do not support chain",
1190         "identifiers. Therefore it is useful to enter a [TT].pdb[tt] file name at",
1191         "the [TT]-o[tt] option when you want to convert a multi-chain [TT].pdb[tt] file.",
1192         "[PAR]",
1193
1194         "The option [TT]-vsite[tt] removes hydrogen and fast improper dihedral",
1195         "motions. Angular and out-of-plane motions can be removed by changing",
1196         "hydrogens into virtual sites and fixing angles, which fixes their",
1197         "position relative to neighboring atoms. Additionally, all atoms in the",
1198         "aromatic rings of the standard amino acids (i.e. PHE, TRP, TYR and HIS)",
1199         "can be converted into virtual sites, eliminating the fast improper dihedral",
1200         "fluctuations in these rings. [BB]Note[bb] that in this case all other hydrogen",
1201         "atoms are also converted to virtual sites. The mass of all atoms that are",
1202         "converted into virtual sites, is added to the heavy atoms.[PAR]",
1203         "Also slowing down of dihedral motion can be done with [TT]-heavyh[tt]",
1204         "done by increasing the hydrogen-mass by a factor of 4. This is also",
1205         "done for water hydrogens to slow down the rotational motion of water.",
1206         "The increase in mass of the hydrogens is subtracted from the bonded",
1207         "(heavy) atom so that the total mass of the system remains the same."
1208     };
1209
1210
1211     FILE             *fp, *top_file, *top_file2, *itp_file = NULL;
1212     int               natom, nres;
1213     t_atoms           pdba_all, *pdba;
1214     t_atoms          *atoms;
1215     t_resinfo        *ri;
1216     t_blocka         *block;
1217     int               chain, nch, maxch, nwaterchain;
1218     t_pdbchain       *pdb_ch;
1219     t_chain          *chains, *cc;
1220     char              select[STRLEN];
1221     int               nincl, nmol;
1222     char            **incls;
1223     t_mols           *mols;
1224     char            **gnames;
1225     int               ePBC;
1226     matrix            box;
1227     rvec              box_space;
1228     int               i, j, k, l, nrtp;
1229     int              *swap_index, si;
1230     t_restp          *restp;
1231     t_hackblock      *ah;
1232     t_symtab          symtab;
1233     gpp_atomtype_t    atype;
1234     gmx_residuetype_t rt;
1235     const char       *top_fn;
1236     char              fn[256], itp_fn[STRLEN], posre_fn[STRLEN], buf_fn[STRLEN];
1237     char              molname[STRLEN], title[STRLEN], quote[STRLEN];
1238     char             *c, forcefield[STRLEN], ffdir[STRLEN];
1239     char              ffname[STRLEN], suffix[STRLEN], buf[STRLEN];
1240     char             *watermodel;
1241     const char       *watres;
1242     int               nrtpf;
1243     char            **rtpf;
1244     char              rtp[STRLEN];
1245     int               nrrn;
1246     char            **rrn;
1247     int               nrtprename, naa;
1248     rtprename_t      *rtprename = NULL;
1249     int               nah, nNtdb, nCtdb, ntdblist;
1250     t_hackblock      *ntdb, *ctdb, **tdblist;
1251     int               nssbonds;
1252     t_ssbond         *ssbonds;
1253     rvec             *pdbx, *x;
1254     gmx_bool          bVsites = FALSE, bWat, bPrevWat = FALSE, bITP, bVsiteAromatics = FALSE, bCheckMerge;
1255     real              mHmult  = 0;
1256     t_hackblock      *hb_chain;
1257     t_restp          *restp_chain;
1258     output_env_t      oenv;
1259     const char       *p_restype;
1260     int               rc;
1261     int               this_atomnum;
1262     int               prev_atomnum;
1263     const char     *  prev_atomname;
1264     const char     *  this_atomname;
1265     const char     *  prev_resname;
1266     const char     *  this_resname;
1267     int               prev_resnum;
1268     int               this_resnum;
1269     char              prev_chainid;
1270     char              this_chainid;
1271     int               prev_chainnumber;
1272     int               this_chainnumber;
1273     int               nid_used;
1274     int               this_chainstart;
1275     int               prev_chainstart;
1276     gmx_bool          bMerged;
1277     int               nchainmerges;
1278
1279     gmx_atomprop_t    aps;
1280
1281     t_filenm          fnm[] = {
1282         { efSTX, "-f", "eiwit.pdb", ffREAD  },
1283         { efSTO, "-o", "conf",      ffWRITE },
1284         { efTOP, NULL, NULL,        ffWRITE },
1285         { efITP, "-i", "posre",     ffWRITE },
1286         { efNDX, "-n", "clean",     ffOPTWR },
1287         { efSTO, "-q", "clean.pdb", ffOPTWR }
1288     };
1289 #define NFILE asize(fnm)
1290
1291
1292     /* Command line arguments must be static */
1293     static gmx_bool    bNewRTP        = FALSE;
1294     static gmx_bool    bInter         = FALSE, bCysMan = FALSE;
1295     static gmx_bool    bLysMan        = FALSE, bAspMan = FALSE, bGluMan = FALSE, bHisMan = FALSE;
1296     static gmx_bool    bGlnMan        = FALSE, bArgMan = FALSE;
1297     static gmx_bool    bTerMan        = FALSE, bUnA = FALSE, bHeavyH;
1298     static gmx_bool    bSort          = TRUE, bAllowMissing = FALSE, bRemoveH = FALSE;
1299     static gmx_bool    bDeuterate     = FALSE, bVerbose = FALSE, bChargeGroups = TRUE, bCmap = TRUE;
1300     static gmx_bool    bRenumRes      = FALSE, bRTPresname = FALSE;
1301     static real        angle          = 135.0, distance = 0.3, posre_fc = 1000;
1302     static real        long_bond_dist = 0.25, short_bond_dist = 0.05;
1303     static const char *vsitestr[]     = { NULL, "none", "hydrogens", "aromatics", NULL };
1304     static const char *watstr[]       = { NULL, "select", "none", "spc", "spce", "tip3p", "tip4p", "tip5p", NULL };
1305     static const char *chainsep[]     = { NULL, "id_or_ter", "id_and_ter", "ter", "id", "interactive", NULL };
1306     static const char *merge[]        = {NULL, "no", "all", "interactive", NULL };
1307     static const char *ff             = "select";
1308
1309     t_pargs            pa[] = {
1310         { "-newrtp", FALSE, etBOOL, {&bNewRTP},
1311           "HIDDENWrite the residue database in new format to [TT]new.rtp[tt]"},
1312         { "-lb",     FALSE, etREAL, {&long_bond_dist},
1313           "HIDDENLong bond warning distance" },
1314         { "-sb",     FALSE, etREAL, {&short_bond_dist},
1315           "HIDDENShort bond warning distance" },
1316         { "-chainsep", FALSE, etENUM, {chainsep},
1317           "Condition in PDB files when a new chain should be started (adding termini)" },
1318         { "-merge",  FALSE, etENUM, {&merge},
1319           "Merge multiple chains into a single [moleculetype]" },
1320         { "-ff",     FALSE, etSTR,  {&ff},
1321           "Force field, interactive by default. Use [TT]-h[tt] for information." },
1322         { "-water",  FALSE, etENUM, {watstr},
1323           "Water model to use" },
1324         { "-inter",  FALSE, etBOOL, {&bInter},
1325           "Set the next 8 options to interactive"},
1326         { "-ss",     FALSE, etBOOL, {&bCysMan},
1327           "Interactive SS bridge selection" },
1328         { "-ter",    FALSE, etBOOL, {&bTerMan},
1329           "Interactive termini selection, instead of charged (default)" },
1330         { "-lys",    FALSE, etBOOL, {&bLysMan},
1331           "Interactive lysine selection, instead of charged" },
1332         { "-arg",    FALSE, etBOOL, {&bArgMan},
1333           "Interactive arginine selection, instead of charged" },
1334         { "-asp",    FALSE, etBOOL, {&bAspMan},
1335           "Interactive aspartic acid selection, instead of charged" },
1336         { "-glu",    FALSE, etBOOL, {&bGluMan},
1337           "Interactive glutamic acid selection, instead of charged" },
1338         { "-gln",    FALSE, etBOOL, {&bGlnMan},
1339           "Interactive glutamine selection, instead of neutral" },
1340         { "-his",    FALSE, etBOOL, {&bHisMan},
1341           "Interactive histidine selection, instead of checking H-bonds" },
1342         { "-angle",  FALSE, etREAL, {&angle},
1343           "Minimum hydrogen-donor-acceptor angle for a H-bond (degrees)" },
1344         { "-dist",   FALSE, etREAL, {&distance},
1345           "Maximum donor-acceptor distance for a H-bond (nm)" },
1346         { "-una",    FALSE, etBOOL, {&bUnA},
1347           "Select aromatic rings with united CH atoms on phenylalanine, "
1348           "tryptophane and tyrosine" },
1349         { "-sort",   FALSE, etBOOL, {&bSort},
1350           "HIDDENSort the residues according to database, turning this off is dangerous as charge groups might be broken in parts" },
1351         { "-ignh",   FALSE, etBOOL, {&bRemoveH},
1352           "Ignore hydrogen atoms that are in the coordinate file" },
1353         { "-missing", FALSE, etBOOL, {&bAllowMissing},
1354           "Continue when atoms are missing, dangerous" },
1355         { "-v",      FALSE, etBOOL, {&bVerbose},
1356           "Be slightly more verbose in messages" },
1357         { "-posrefc", FALSE, etREAL, {&posre_fc},
1358           "Force constant for position restraints" },
1359         { "-vsite",  FALSE, etENUM, {vsitestr},
1360           "Convert atoms to virtual sites" },
1361         { "-heavyh", FALSE, etBOOL, {&bHeavyH},
1362           "Make hydrogen atoms heavy" },
1363         { "-deuterate", FALSE, etBOOL, {&bDeuterate},
1364           "Change the mass of hydrogens to 2 amu" },
1365         { "-chargegrp", TRUE, etBOOL, {&bChargeGroups},
1366           "Use charge groups in the [TT].rtp[tt] file"  },
1367         { "-cmap", TRUE, etBOOL, {&bCmap},
1368           "Use cmap torsions (if enabled in the [TT].rtp[tt] file)"  },
1369         { "-renum", TRUE, etBOOL, {&bRenumRes},
1370           "Renumber the residues consecutively in the output"  },
1371         { "-rtpres", TRUE, etBOOL, {&bRTPresname},
1372           "Use [TT].rtp[tt] entry names as residue names"  }
1373     };
1374 #define NPARGS asize(pa)
1375
1376     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc,
1377                            0, NULL, &oenv))
1378     {
1379         return 0;
1380     }
1381
1382     /* Force field selection, interactive or direct */
1383     choose_ff(strcmp(ff, "select") == 0 ? NULL : ff,
1384               forcefield, sizeof(forcefield),
1385               ffdir, sizeof(ffdir));
1386
1387     if (strlen(forcefield) > 0)
1388     {
1389         strcpy(ffname, forcefield);
1390         ffname[0] = toupper(ffname[0]);
1391     }
1392     else
1393     {
1394         gmx_fatal(FARGS, "Empty forcefield string");
1395     }
1396
1397     printf("\nUsing the %s force field in directory %s\n\n",
1398            ffname, ffdir);
1399
1400     choose_watermodel(watstr[0], ffdir, &watermodel);
1401
1402     if (bInter)
1403     {
1404         /* if anything changes here, also change description of -inter */
1405         bCysMan = TRUE;
1406         bTerMan = TRUE;
1407         bLysMan = TRUE;
1408         bArgMan = TRUE;
1409         bAspMan = TRUE;
1410         bGluMan = TRUE;
1411         bGlnMan = TRUE;
1412         bHisMan = TRUE;
1413     }
1414
1415     if (bHeavyH)
1416     {
1417         mHmult = 4.0;
1418     }
1419     else if (bDeuterate)
1420     {
1421         mHmult = 2.0;
1422     }
1423     else
1424     {
1425         mHmult = 1.0;
1426     }
1427
1428     switch (vsitestr[0][0])
1429     {
1430         case 'n': /* none */
1431             bVsites         = FALSE;
1432             bVsiteAromatics = FALSE;
1433             break;
1434         case 'h': /* hydrogens */
1435             bVsites         = TRUE;
1436             bVsiteAromatics = FALSE;
1437             break;
1438         case 'a': /* aromatics */
1439             bVsites         = TRUE;
1440             bVsiteAromatics = TRUE;
1441             break;
1442         default:
1443             gmx_fatal(FARGS, "DEATH HORROR in $s (%d): vsitestr[0]='%s'",
1444                       __FILE__, __LINE__, vsitestr[0]);
1445     } /* end switch */
1446
1447     /* Open the symbol table */
1448     open_symtab(&symtab);
1449
1450     /* Residue type database */
1451     gmx_residuetype_init(&rt);
1452
1453     /* Read residue renaming database(s), if present */
1454     nrrn = fflib_search_file_end(ffdir, ".r2b", FALSE, &rrn);
1455
1456     nrtprename = 0;
1457     rtprename  = NULL;
1458     for (i = 0; i < nrrn; i++)
1459     {
1460         fp = fflib_open(rrn[i]);
1461         read_rtprename(rrn[i], fp, &nrtprename, &rtprename);
1462         gmx_ffclose(fp);
1463         sfree(rrn[i]);
1464     }
1465     sfree(rrn);
1466
1467     /* Add all alternative names from the residue renaming database to the list of recognized amino/nucleic acids. */
1468     naa = 0;
1469     for (i = 0; i < nrtprename; i++)
1470     {
1471         rc = gmx_residuetype_get_type(rt, rtprename[i].gmx, &p_restype);
1472
1473         /* Only add names if the 'standard' gromacs/iupac base name was found */
1474         if (rc == 0)
1475         {
1476             gmx_residuetype_add(rt, rtprename[i].main, p_restype);
1477             gmx_residuetype_add(rt, rtprename[i].nter, p_restype);
1478             gmx_residuetype_add(rt, rtprename[i].cter, p_restype);
1479             gmx_residuetype_add(rt, rtprename[i].bter, p_restype);
1480         }
1481     }
1482
1483     clear_mat(box);
1484     if (watermodel != NULL && (strstr(watermodel, "4p") ||
1485                                strstr(watermodel, "4P")))
1486     {
1487         watres = "HO4";
1488     }
1489     else if (watermodel != NULL && (strstr(watermodel, "5p") ||
1490                                     strstr(watermodel, "5P")))
1491     {
1492         watres = "HO5";
1493     }
1494     else
1495     {
1496         watres = "HOH";
1497     }
1498
1499     aps   = gmx_atomprop_init();
1500     natom = read_pdball(opt2fn("-f", NFILE, fnm), opt2fn_null("-q", NFILE, fnm), title,
1501                         &pdba_all, &pdbx, &ePBC, box, bRemoveH, &symtab, rt, watres,
1502                         aps, bVerbose);
1503
1504     if (natom == 0)
1505     {
1506         gmx_fatal(FARGS, "No atoms found in pdb file %s\n", opt2fn("-f", NFILE, fnm));
1507     }
1508
1509     printf("Analyzing pdb file\n");
1510     nch         = 0;
1511     maxch       = 0;
1512     nwaterchain = 0;
1513
1514     modify_chain_numbers(&pdba_all, chainsep[0]);
1515
1516     nchainmerges        = 0;
1517
1518     this_atomname       = NULL;
1519     this_atomnum        = -1;
1520     this_resname        = NULL;
1521     this_resnum         = -1;
1522     this_chainid        = '?';
1523     this_chainnumber    = -1;
1524     this_chainstart     = 0;
1525     /* Keep the compiler happy */
1526     prev_chainstart     = 0;
1527
1528     pdb_ch = NULL;
1529
1530     bMerged = FALSE;
1531     for (i = 0; (i < natom); i++)
1532     {
1533         ri = &pdba_all.resinfo[pdba_all.atom[i].resind];
1534
1535         prev_atomname      = this_atomname;
1536         prev_atomnum       = this_atomnum;
1537         prev_resname       = this_resname;
1538         prev_resnum        = this_resnum;
1539         prev_chainid       = this_chainid;
1540         prev_chainnumber   = this_chainnumber;
1541         if (!bMerged)
1542         {
1543             prev_chainstart    = this_chainstart;
1544         }
1545
1546         this_atomname      = *pdba_all.atomname[i];
1547         this_atomnum       = (pdba_all.pdbinfo != NULL) ? pdba_all.pdbinfo[i].atomnr : i+1;
1548         this_resname       = *ri->name;
1549         this_resnum        = ri->nr;
1550         this_chainid       = ri->chainid;
1551         this_chainnumber   = ri->chainnum;
1552
1553         bWat = gmx_strcasecmp(*ri->name, watres) == 0;
1554         if ((i == 0) || (this_chainnumber != prev_chainnumber) || (bWat != bPrevWat))
1555         {
1556             this_chainstart = pdba_all.atom[i].resind;
1557
1558             bMerged = FALSE;
1559             if (i > 0 && !bWat)
1560             {
1561                 if (!strncmp(merge[0], "int", 3))
1562                 {
1563                     printf("Merge chain ending with residue %s%d (chain id '%c', atom %d %s) and chain starting with\n"
1564                            "residue %s%d (chain id '%c', atom %d %s) into a single moleculetype (keeping termini)? [n/y]\n",
1565                            prev_resname, prev_resnum, prev_chainid, prev_atomnum, prev_atomname,
1566                            this_resname, this_resnum, this_chainid, this_atomnum, this_atomname);
1567
1568                     if (NULL == fgets(select, STRLEN-1, stdin))
1569                     {
1570                         gmx_fatal(FARGS, "Error reading from stdin");
1571                     }
1572                     bMerged = (select[0] == 'y');
1573                 }
1574                 else if (!strncmp(merge[0], "all", 3))
1575                 {
1576                     bMerged = TRUE;
1577                 }
1578             }
1579
1580             if (bMerged)
1581             {
1582                 pdb_ch[nch-1].chainstart[pdb_ch[nch-1].nterpairs] =
1583                     pdba_all.atom[i].resind - prev_chainstart;
1584                 pdb_ch[nch-1].nterpairs++;
1585                 srenew(pdb_ch[nch-1].chainstart, pdb_ch[nch-1].nterpairs+1);
1586                 nchainmerges++;
1587             }
1588             else
1589             {
1590                 /* set natom for previous chain */
1591                 if (nch > 0)
1592                 {
1593                     pdb_ch[nch-1].natom = i-pdb_ch[nch-1].start;
1594                 }
1595                 if (bWat)
1596                 {
1597                     nwaterchain++;
1598                     ri->chainid = ' ';
1599                 }
1600                 /* check if chain identifier was used before */
1601                 for (j = 0; (j < nch); j++)
1602                 {
1603                     if (pdb_ch[j].chainid != ' ' && pdb_ch[j].chainid == ri->chainid)
1604                     {
1605                         printf("WARNING: Chain identifier '%c' is used in two non-sequential blocks.\n"
1606                                "They will be treated as separate chains unless you reorder your file.\n",
1607                                ri->chainid);
1608                     }
1609                 }
1610                 if (nch == maxch)
1611                 {
1612                     maxch += 16;
1613                     srenew(pdb_ch, maxch);
1614                 }
1615                 pdb_ch[nch].chainid  = ri->chainid;
1616                 pdb_ch[nch].chainnum = ri->chainnum;
1617                 pdb_ch[nch].start    = i;
1618                 pdb_ch[nch].bAllWat  = bWat;
1619                 if (bWat)
1620                 {
1621                     pdb_ch[nch].nterpairs = 0;
1622                 }
1623                 else
1624                 {
1625                     pdb_ch[nch].nterpairs = 1;
1626                 }
1627                 snew(pdb_ch[nch].chainstart, pdb_ch[nch].nterpairs+1);
1628                 /* modified [nch] to [0] below */
1629                 pdb_ch[nch].chainstart[0] = 0;
1630                 nch++;
1631             }
1632         }
1633         bPrevWat = bWat;
1634     }
1635     pdb_ch[nch-1].natom = natom-pdb_ch[nch-1].start;
1636
1637     /* set all the water blocks at the end of the chain */
1638     snew(swap_index, nch);
1639     j = 0;
1640     for (i = 0; i < nch; i++)
1641     {
1642         if (!pdb_ch[i].bAllWat)
1643         {
1644             swap_index[j] = i;
1645             j++;
1646         }
1647     }
1648     for (i = 0; i < nch; i++)
1649     {
1650         if (pdb_ch[i].bAllWat)
1651         {
1652             swap_index[j] = i;
1653             j++;
1654         }
1655     }
1656     if (nwaterchain > 1)
1657     {
1658         printf("Moved all the water blocks to the end\n");
1659     }
1660
1661     snew(chains, nch);
1662     /* copy pdb data and x for all chains */
1663     for (i = 0; (i < nch); i++)
1664     {
1665         si                   = swap_index[i];
1666         chains[i].chainid    = pdb_ch[si].chainid;
1667         chains[i].chainnum   = pdb_ch[si].chainnum;
1668         chains[i].bAllWat    = pdb_ch[si].bAllWat;
1669         chains[i].nterpairs  = pdb_ch[si].nterpairs;
1670         chains[i].chainstart = pdb_ch[si].chainstart;
1671         snew(chains[i].ntdb, pdb_ch[si].nterpairs);
1672         snew(chains[i].ctdb, pdb_ch[si].nterpairs);
1673         snew(chains[i].r_start, pdb_ch[si].nterpairs);
1674         snew(chains[i].r_end, pdb_ch[si].nterpairs);
1675
1676         snew(chains[i].pdba, 1);
1677         init_t_atoms(chains[i].pdba, pdb_ch[si].natom, TRUE);
1678         snew(chains[i].x, chains[i].pdba->nr);
1679         for (j = 0; j < chains[i].pdba->nr; j++)
1680         {
1681             chains[i].pdba->atom[j] = pdba_all.atom[pdb_ch[si].start+j];
1682             snew(chains[i].pdba->atomname[j], 1);
1683             *chains[i].pdba->atomname[j] =
1684                 strdup(*pdba_all.atomname[pdb_ch[si].start+j]);
1685             chains[i].pdba->pdbinfo[j] = pdba_all.pdbinfo[pdb_ch[si].start+j];
1686             copy_rvec(pdbx[pdb_ch[si].start+j], chains[i].x[j]);
1687         }
1688         /* Re-index the residues assuming that the indices are continuous */
1689         k                    = chains[i].pdba->atom[0].resind;
1690         nres                 = chains[i].pdba->atom[chains[i].pdba->nr-1].resind - k + 1;
1691         chains[i].pdba->nres = nres;
1692         for (j = 0; j < chains[i].pdba->nr; j++)
1693         {
1694             chains[i].pdba->atom[j].resind -= k;
1695         }
1696         srenew(chains[i].pdba->resinfo, nres);
1697         for (j = 0; j < nres; j++)
1698         {
1699             chains[i].pdba->resinfo[j] = pdba_all.resinfo[k+j];
1700             snew(chains[i].pdba->resinfo[j].name, 1);
1701             *chains[i].pdba->resinfo[j].name = strdup(*pdba_all.resinfo[k+j].name);
1702             /* make all chain identifiers equal to that of the chain */
1703             chains[i].pdba->resinfo[j].chainid = pdb_ch[si].chainid;
1704         }
1705     }
1706
1707     if (nchainmerges > 0)
1708     {
1709         printf("\nMerged chains into joint molecule definitions at %d places.\n\n",
1710                nchainmerges);
1711     }
1712
1713     printf("There are %d chains and %d blocks of water and "
1714            "%d residues with %d atoms\n",
1715            nch-nwaterchain, nwaterchain,
1716            pdba_all.resinfo[pdba_all.atom[natom-1].resind].nr, natom);
1717
1718     printf("\n  %5s  %4s %6s\n", "chain", "#res", "#atoms");
1719     for (i = 0; (i < nch); i++)
1720     {
1721         printf("  %d '%c' %5d %6d  %s\n",
1722                i+1, chains[i].chainid ? chains[i].chainid : '-',
1723                chains[i].pdba->nres, chains[i].pdba->nr,
1724                chains[i].bAllWat ? "(only water)" : "");
1725     }
1726     printf("\n");
1727
1728     check_occupancy(&pdba_all, opt2fn("-f", NFILE, fnm), bVerbose);
1729
1730     /* Read atomtypes... */
1731     atype = read_atype(ffdir, &symtab);
1732
1733     /* read residue database */
1734     printf("Reading residue database... (%s)\n", forcefield);
1735     nrtpf = fflib_search_file_end(ffdir, ".rtp", TRUE, &rtpf);
1736     nrtp  = 0;
1737     restp = NULL;
1738     for (i = 0; i < nrtpf; i++)
1739     {
1740         read_resall(rtpf[i], &nrtp, &restp, atype, &symtab, FALSE);
1741         sfree(rtpf[i]);
1742     }
1743     sfree(rtpf);
1744     if (bNewRTP)
1745     {
1746         /* Not correct with multiple rtp input files with different bonded types */
1747         fp = gmx_fio_fopen("new.rtp", "w");
1748         print_resall(fp, nrtp, restp, atype);
1749         gmx_fio_fclose(fp);
1750     }
1751
1752     /* read hydrogen database */
1753     nah = read_h_db(ffdir, &ah);
1754
1755     /* Read Termini database... */
1756     nNtdb = read_ter_db(ffdir, 'n', &ntdb, atype);
1757     nCtdb = read_ter_db(ffdir, 'c', &ctdb, atype);
1758
1759     top_fn   = ftp2fn(efTOP, NFILE, fnm);
1760     top_file = gmx_fio_fopen(top_fn, "w");
1761
1762     print_top_header(top_file, top_fn, FALSE, ffdir, mHmult);
1763
1764     nincl = 0;
1765     nmol  = 0;
1766     incls = NULL;
1767     mols  = NULL;
1768     nres  = 0;
1769     for (chain = 0; (chain < nch); chain++)
1770     {
1771         cc = &(chains[chain]);
1772
1773         /* set pdba, natom and nres to the current chain */
1774         pdba  = cc->pdba;
1775         x     = cc->x;
1776         natom = cc->pdba->nr;
1777         nres  = cc->pdba->nres;
1778
1779         if (cc->chainid && ( cc->chainid != ' ' ) )
1780         {
1781             printf("Processing chain %d '%c' (%d atoms, %d residues)\n",
1782                    chain+1, cc->chainid, natom, nres);
1783         }
1784         else
1785         {
1786             printf("Processing chain %d (%d atoms, %d residues)\n",
1787                    chain+1, natom, nres);
1788         }
1789
1790         process_chain(pdba, x, bUnA, bUnA, bUnA, bLysMan, bAspMan, bGluMan,
1791                       bHisMan, bArgMan, bGlnMan, angle, distance, &symtab,
1792                       nrtprename, rtprename);
1793
1794         cc->chainstart[cc->nterpairs] = pdba->nres;
1795         j = 0;
1796         for (i = 0; i < cc->nterpairs; i++)
1797         {
1798             find_nc_ter(pdba, cc->chainstart[i], cc->chainstart[i+1],
1799                         &(cc->r_start[j]), &(cc->r_end[j]), rt);
1800
1801             if (cc->r_start[j] >= 0 && cc->r_end[j] >= 0)
1802             {
1803                 j++;
1804             }
1805         }
1806         cc->nterpairs = j;
1807         if (cc->nterpairs == 0)
1808         {
1809             printf("Problem with chain definition, or missing terminal residues.\n"
1810                    "This chain does not appear to contain a recognized chain molecule.\n"
1811                    "If this is incorrect, you can edit residuetypes.dat to modify the behavior.\n");
1812         }
1813
1814         /* Check for disulfides and other special bonds */
1815         nssbonds = mk_specbonds(pdba, x, bCysMan, &ssbonds, bVerbose);
1816
1817         if (nrtprename > 0)
1818         {
1819             rename_resrtp(pdba, cc->nterpairs, cc->r_start, cc->r_end, nrtprename, rtprename,
1820                           &symtab, bVerbose);
1821         }
1822
1823         if (debug)
1824         {
1825             if (nch == 1)
1826             {
1827                 sprintf(fn, "chain.pdb");
1828             }
1829             else
1830             {
1831                 sprintf(fn, "chain_%c%d.pdb", cc->chainid, cc->chainnum);
1832             }
1833             write_sto_conf(fn, title, pdba, x, NULL, ePBC, box);
1834         }
1835
1836
1837         for (i = 0; i < cc->nterpairs; i++)
1838         {
1839
1840             /* Set termini.
1841              * We first apply a filter so we only have the
1842              * termini that can be applied to the residue in question
1843              * (or a generic terminus if no-residue specific is available).
1844              */
1845             /* First the N terminus */
1846             if (nNtdb > 0)
1847             {
1848                 tdblist = filter_ter(nrtp, restp, nNtdb, ntdb,
1849                                      *pdba->resinfo[cc->r_start[i]].name,
1850                                      *pdba->resinfo[cc->r_start[i]].rtp,
1851                                      &ntdblist);
1852                 if (ntdblist == 0)
1853                 {
1854                     printf("No suitable end (N or 5') terminus found in database - assuming this residue\n"
1855                            "is already in a terminus-specific form and skipping terminus selection.\n");
1856                     cc->ntdb[i] = NULL;
1857                 }
1858                 else
1859                 {
1860                     if (bTerMan && ntdblist > 1)
1861                     {
1862                         sprintf(select, "Select start terminus type for %s-%d",
1863                                 *pdba->resinfo[cc->r_start[i]].name,
1864                                 pdba->resinfo[cc->r_start[i]].nr);
1865                         cc->ntdb[i] = choose_ter(ntdblist, tdblist, select);
1866                     }
1867                     else
1868                     {
1869                         cc->ntdb[i] = tdblist[0];
1870                     }
1871
1872                     printf("Start terminus %s-%d: %s\n",
1873                            *pdba->resinfo[cc->r_start[i]].name,
1874                            pdba->resinfo[cc->r_start[i]].nr,
1875                            (cc->ntdb[i])->name);
1876                     sfree(tdblist);
1877                 }
1878             }
1879             else
1880             {
1881                 cc->ntdb[i] = NULL;
1882             }
1883
1884             /* And the C terminus */
1885             if (nCtdb > 0)
1886             {
1887                 tdblist = filter_ter(nrtp, restp, nCtdb, ctdb,
1888                                      *pdba->resinfo[cc->r_end[i]].name,
1889                                      *pdba->resinfo[cc->r_end[i]].rtp,
1890                                      &ntdblist);
1891                 if (ntdblist == 0)
1892                 {
1893                     printf("No suitable end (C or 3') terminus found in database - assuming this residue\n"
1894                            "is already in a terminus-specific form and skipping terminus selection.\n");
1895                     cc->ctdb[i] = NULL;
1896                 }
1897                 else
1898                 {
1899                     if (bTerMan && ntdblist > 1)
1900                     {
1901                         sprintf(select, "Select end terminus type for %s-%d",
1902                                 *pdba->resinfo[cc->r_end[i]].name,
1903                                 pdba->resinfo[cc->r_end[i]].nr);
1904                         cc->ctdb[i] = choose_ter(ntdblist, tdblist, select);
1905                     }
1906                     else
1907                     {
1908                         cc->ctdb[i] = tdblist[0];
1909                     }
1910                     printf("End terminus %s-%d: %s\n",
1911                            *pdba->resinfo[cc->r_end[i]].name,
1912                            pdba->resinfo[cc->r_end[i]].nr,
1913                            (cc->ctdb[i])->name);
1914                     sfree(tdblist);
1915                 }
1916             }
1917             else
1918             {
1919                 cc->ctdb[i] = NULL;
1920             }
1921         }
1922         /* lookup hackblocks and rtp for all residues */
1923         get_hackblocks_rtp(&hb_chain, &restp_chain,
1924                            nrtp, restp, pdba->nres, pdba->resinfo,
1925                            cc->nterpairs, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end);
1926         /* ideally, now we would not need the rtp itself anymore, but do
1927            everything using the hb and restp arrays. Unfortunately, that
1928            requires some re-thinking of code in gen_vsite.c, which I won't
1929            do now :( AF 26-7-99 */
1930
1931         rename_atoms(NULL, ffdir,
1932                      pdba, &symtab, restp_chain, FALSE, rt, FALSE, bVerbose);
1933
1934         match_atomnames_with_rtp(restp_chain, hb_chain, pdba, x, bVerbose);
1935
1936         if (bSort)
1937         {
1938             block = new_blocka();
1939             snew(gnames, 1);
1940             sort_pdbatoms(restp_chain, natom, &pdba, &x, block, &gnames);
1941             natom = remove_duplicate_atoms(pdba, x, bVerbose);
1942             if (ftp2bSet(efNDX, NFILE, fnm))
1943             {
1944                 if (bRemoveH)
1945                 {
1946                     fprintf(stderr, "WARNING: with the -remh option the generated "
1947                             "index file (%s) might be useless\n"
1948                             "(the index file is generated before hydrogens are added)",
1949                             ftp2fn(efNDX, NFILE, fnm));
1950                 }
1951                 write_index(ftp2fn(efNDX, NFILE, fnm), block, gnames, FALSE, 0);
1952             }
1953             for (i = 0; i < block->nr; i++)
1954             {
1955                 sfree(gnames[i]);
1956             }
1957             sfree(gnames);
1958             done_blocka(block);
1959         }
1960         else
1961         {
1962             fprintf(stderr, "WARNING: "
1963                     "without sorting no check for duplicate atoms can be done\n");
1964         }
1965
1966         /* Generate Hydrogen atoms (and termini) in the sequence */
1967         printf("Generating any missing hydrogen atoms and/or adding termini.\n");
1968         natom = add_h(&pdba, &x, nah, ah,
1969                       cc->nterpairs, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end, bAllowMissing,
1970                       NULL, NULL, TRUE, FALSE);
1971         printf("Now there are %d residues with %d atoms\n",
1972                pdba->nres, pdba->nr);
1973         if (debug)
1974         {
1975             write_pdbfile(debug, title, pdba, x, ePBC, box, ' ', 0, NULL, TRUE);
1976         }
1977
1978         if (debug)
1979         {
1980             for (i = 0; (i < natom); i++)
1981             {
1982                 fprintf(debug, "Res %s%d atom %d %s\n",
1983                         *(pdba->resinfo[pdba->atom[i].resind].name),
1984                         pdba->resinfo[pdba->atom[i].resind].nr, i+1, *pdba->atomname[i]);
1985             }
1986         }
1987
1988         strcpy(posre_fn, ftp2fn(efITP, NFILE, fnm));
1989
1990         /* make up molecule name(s) */
1991
1992         k = (cc->nterpairs > 0 && cc->r_start[0] >= 0) ? cc->r_start[0] : 0;
1993
1994         gmx_residuetype_get_type(rt, *pdba->resinfo[k].name, &p_restype);
1995
1996         suffix[0] = '\0';
1997
1998         if (cc->bAllWat)
1999         {
2000             sprintf(molname, "Water");
2001         }
2002         else
2003         {
2004             this_chainid = cc->chainid;
2005
2006             /* Add the chain id if we have one */
2007             if (this_chainid != ' ')
2008             {
2009                 sprintf(buf, "_chain_%c", this_chainid);
2010                 strcat(suffix, buf);
2011             }
2012
2013             /* Check if there have been previous chains with the same id */
2014             nid_used = 0;
2015             for (k = 0; k < chain; k++)
2016             {
2017                 if (cc->chainid == chains[k].chainid)
2018                 {
2019                     nid_used++;
2020                 }
2021             }
2022             /* Add the number for this chain identifier if there are multiple copies */
2023             if (nid_used > 0)
2024             {
2025
2026                 sprintf(buf, "%d", nid_used+1);
2027                 strcat(suffix, buf);
2028             }
2029
2030             if (strlen(suffix) > 0)
2031             {
2032                 sprintf(molname, "%s%s", p_restype, suffix);
2033             }
2034             else
2035             {
2036                 strcpy(molname, p_restype);
2037             }
2038         }
2039
2040         if ((nch-nwaterchain > 1) && !cc->bAllWat)
2041         {
2042             bITP = TRUE;
2043             strcpy(itp_fn, top_fn);
2044             printf("Chain time...\n");
2045             c = strrchr(itp_fn, '.');
2046             sprintf(c, "_%s.itp", molname);
2047             c = strrchr(posre_fn, '.');
2048             sprintf(c, "_%s.itp", molname);
2049             if (strcmp(itp_fn, posre_fn) == 0)
2050             {
2051                 strcpy(buf_fn, posre_fn);
2052                 c  = strrchr(buf_fn, '.');
2053                 *c = '\0';
2054                 sprintf(posre_fn, "%s_pr.itp", buf_fn);
2055             }
2056
2057             nincl++;
2058             srenew(incls, nincl);
2059             incls[nincl-1] = strdup(itp_fn);
2060             itp_file       = gmx_fio_fopen(itp_fn, "w");
2061         }
2062         else
2063         {
2064             bITP = FALSE;
2065         }
2066
2067         srenew(mols, nmol+1);
2068         if (cc->bAllWat)
2069         {
2070             mols[nmol].name = strdup("SOL");
2071             mols[nmol].nr   = pdba->nres;
2072         }
2073         else
2074         {
2075             mols[nmol].name = strdup(molname);
2076             mols[nmol].nr   = 1;
2077         }
2078         nmol++;
2079
2080         if (bITP)
2081         {
2082             print_top_comment(itp_file, itp_fn, ffdir, TRUE);
2083         }
2084
2085         if (cc->bAllWat)
2086         {
2087             top_file2 = NULL;
2088         }
2089         else
2090         if (bITP)
2091         {
2092             top_file2 = itp_file;
2093         }
2094         else
2095         {
2096             top_file2 = top_file;
2097         }
2098
2099         pdb2top(top_file2, posre_fn, molname, pdba, &x, atype, &symtab,
2100                 nrtp, restp,
2101                 restp_chain, hb_chain,
2102                 bAllowMissing,
2103                 bVsites, bVsiteAromatics, ffdir,
2104                 mHmult, nssbonds, ssbonds,
2105                 long_bond_dist, short_bond_dist, bDeuterate, bChargeGroups, bCmap,
2106                 bRenumRes, bRTPresname);
2107
2108         if (!cc->bAllWat)
2109         {
2110             write_posres(posre_fn, pdba, posre_fc);
2111         }
2112
2113         if (bITP)
2114         {
2115             gmx_fio_fclose(itp_file);
2116         }
2117
2118         /* pdba and natom have been reassigned somewhere so: */
2119         cc->pdba = pdba;
2120         cc->x    = x;
2121
2122         if (debug)
2123         {
2124             if (cc->chainid == ' ')
2125             {
2126                 sprintf(fn, "chain.pdb");
2127             }
2128             else
2129             {
2130                 sprintf(fn, "chain_%c.pdb", cc->chainid);
2131             }
2132             cool_quote(quote, 255, NULL);
2133             write_sto_conf(fn, quote, pdba, x, NULL, ePBC, box);
2134         }
2135     }
2136
2137     if (watermodel == NULL)
2138     {
2139         for (chain = 0; chain < nch; chain++)
2140         {
2141             if (chains[chain].bAllWat)
2142             {
2143                 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.");
2144             }
2145         }
2146     }
2147     else
2148     {
2149         sprintf(buf_fn, "%s%c%s.itp", ffdir, DIR_SEPARATOR, watermodel);
2150         if (!fflib_fexist(buf_fn))
2151         {
2152             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.",
2153                       buf_fn, watermodel);
2154         }
2155     }
2156
2157     print_top_mols(top_file, title, ffdir, watermodel, nincl, incls, nmol, mols);
2158     gmx_fio_fclose(top_file);
2159
2160     gmx_residuetype_destroy(rt);
2161
2162     /* now merge all chains back together */
2163     natom = 0;
2164     nres  = 0;
2165     for (i = 0; (i < nch); i++)
2166     {
2167         natom += chains[i].pdba->nr;
2168         nres  += chains[i].pdba->nres;
2169     }
2170     snew(atoms, 1);
2171     init_t_atoms(atoms, natom, FALSE);
2172     for (i = 0; i < atoms->nres; i++)
2173     {
2174         sfree(atoms->resinfo[i].name);
2175     }
2176     sfree(atoms->resinfo);
2177     atoms->nres = nres;
2178     snew(atoms->resinfo, nres);
2179     snew(x, natom);
2180     k = 0;
2181     l = 0;
2182     for (i = 0; (i < nch); i++)
2183     {
2184         if (nch > 1)
2185         {
2186             printf("Including chain %d in system: %d atoms %d residues\n",
2187                    i+1, chains[i].pdba->nr, chains[i].pdba->nres);
2188         }
2189         for (j = 0; (j < chains[i].pdba->nr); j++)
2190         {
2191             atoms->atom[k]         = chains[i].pdba->atom[j];
2192             atoms->atom[k].resind += l; /* l is processed nr of residues */
2193             atoms->atomname[k]     = chains[i].pdba->atomname[j];
2194             atoms->resinfo[atoms->atom[k].resind].chainid = chains[i].chainid;
2195             copy_rvec(chains[i].x[j], x[k]);
2196             k++;
2197         }
2198         for (j = 0; (j < chains[i].pdba->nres); j++)
2199         {
2200             atoms->resinfo[l] = chains[i].pdba->resinfo[j];
2201             if (bRTPresname)
2202             {
2203                 atoms->resinfo[l].name = atoms->resinfo[l].rtp;
2204             }
2205             l++;
2206         }
2207     }
2208
2209     if (nch > 1)
2210     {
2211         fprintf(stderr, "Now there are %d atoms and %d residues\n", k, l);
2212         print_sums(atoms, TRUE);
2213     }
2214
2215     fprintf(stderr, "\nWriting coordinate file...\n");
2216     clear_rvec(box_space);
2217     if (box[0][0] == 0)
2218     {
2219         make_new_box(atoms->nr, x, box, box_space, FALSE);
2220     }
2221     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), title, atoms, x, NULL, ePBC, box);
2222
2223     printf("\t\t--------- PLEASE NOTE ------------\n");
2224     printf("You have successfully generated a topology from: %s.\n",
2225            opt2fn("-f", NFILE, fnm));
2226     if (watermodel != NULL)
2227     {
2228         printf("The %s force field and the %s water model are used.\n",
2229                ffname, watermodel);
2230     }
2231     else
2232     {
2233         printf("The %s force field is used.\n",
2234                ffname);
2235     }
2236     printf("\t\t--------- ETON ESAELP ------------\n");
2237
2238     return 0;
2239 }