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