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