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