Added documentation about proper use of pdb2gmx -ter.

[gromacs/rigid-bodies.git] / src / kernel / pdb2gmx.c

/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
 *
 *
 *                This source code is part of
 * 
 *                 G   R   O   M   A   C   S
 * 
 *          GROningen MAchine for Chemical Simulations
 * 
 *                        VERSION 3.2.0
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team,
 * check out http://www.gromacs.org for more information.

 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 * 
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 * 
 * For more info, check our website at http://www.gromacs.org
 * 
 * And Hey:
 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <time.h>
#include <ctype.h>
#include "sysstuff.h"
#include "typedefs.h"
#include "gmxfio.h"
#include "smalloc.h"
#include "copyrite.h"
#include "string2.h"
#include "confio.h"
#include "symtab.h"
#include "vec.h"
#include "statutil.h"
#include "futil.h"
#include "gmx_fatal.h"
#include "pdbio.h"
#include "toputil.h"
#include "h_db.h"
#include "physics.h"
#include "pgutil.h"
#include "calch.h"
#include "resall.h"
#include "pdb2top.h"
#include "ter_db.h"
#include "strdb.h"
#include "gbutil.h"
#include "genhydro.h"
#include "readinp.h"
#include "atomprop.h"
#include "xlate.h"
#include "specbond.h"
#include "index.h"
#include "hizzie.h"
#include "fflibutil.h"


typedef struct {
  char gmx[6];
  char main[6];
  char nter[6];
  char cter[6];
  char bter[6];
} rtprename_t;


static const char *res2bb_notermini(const char *name,
                                    int nrr,const rtprename_t *rr)
{
  /* NOTE: This function returns the main building block name,
   *       it does not take terminal renaming into account.
   */
  int i;

  i = 0;
  while (i < nrr && gmx_strcasecmp(name,rr[i].gmx) != 0) {
    i++;
  }

  return (i < nrr ? rr[i].main : name);
}

static const char *select_res(int nr,int resnr,
                              const char *name[],const char *expl[],
                              const char *title,
                              int nrr,const rtprename_t *rr)
{
  int sel=0;

  printf("Which %s type do you want for residue %d\n",title,resnr+1);
  for(sel=0; (sel < nr); sel++) {
    printf("%d. %s (%s)\n",
           sel,expl[sel],res2bb_notermini(name[sel],nrr,rr));
  }
  printf("\nType a number:"); fflush(stdout);

  if (scanf("%d",&sel) != 1)
    gmx_fatal(FARGS,"Answer me for res %s %d!",title,resnr+1);
  
  return name[sel];
}

static const char *get_asptp(int resnr,int nrr,const rtprename_t *rr)
{
  enum { easp, easpH, easpNR };
  const char *lh[easpNR] = { "ASP", "ASPH" };
  const char *expl[easpNR] = {
    "Not protonated (charge -1)",
    "Protonated (charge 0)"
  };

  return select_res(easpNR,resnr,lh,expl,"ASPARTIC ACID",nrr,rr);
}

static const char *get_glutp(int resnr,int nrr,const rtprename_t *rr)
{
  enum { eglu, egluH, egluNR };
  const char *lh[egluNR] = { "GLU", "GLUH" };
  const char *expl[egluNR] = {
    "Not protonated (charge -1)",
    "Protonated (charge 0)"
  };

  return select_res(egluNR,resnr,lh,expl,"GLUTAMIC ACID",nrr,rr);
}

static const char *get_glntp(int resnr,int nrr,const rtprename_t *rr)
{
  enum { egln, eglnH, eglnNR };
  const char *lh[eglnNR] = { "GLN", "QLN" };
  const char *expl[eglnNR] = {
    "Not protonated (charge 0)",
    "Protonated (charge +1)"
  };

  return select_res(eglnNR,resnr,lh,expl,"GLUTAMINE",nrr,rr);
}

static const char *get_lystp(int resnr,int nrr,const rtprename_t *rr)
{
  enum { elys, elysH, elysNR };
  const  char *lh[elysNR] = { "LYSN", "LYS" };
  const char *expl[elysNR] = {
    "Not protonated (charge 0)",
    "Protonated (charge +1)"
  };

  return select_res(elysNR,resnr,lh,expl,"LYSINE",nrr,rr);
}

static const char *get_argtp(int resnr,int nrr,const rtprename_t *rr)
{
  enum { earg, eargH, eargNR };
  const  char *lh[eargNR] = { "ARGN", "ARG" };
  const char *expl[eargNR] = {
    "Not protonated (charge 0)",
    "Protonated (charge +1)"
  };

  return select_res(eargNR,resnr,lh,expl,"ARGININE",nrr,rr);
}

static const char *get_cystp(int resnr,int nrr,const rtprename_t *rr)
{
  enum { ecys, ecysH, ecysNR };
  const char *lh[ecysNR] = { "CYS2", "CYS" };
  const char *expl[ecysNR] = {
    "Cysteine in disulfide bridge",
    "Protonated"
  };

  return select_res(ecysNR,resnr,lh,expl,"CYSTEINE",nrr,rr);

}

static const char *get_histp(int resnr,int nrr,const rtprename_t *rr)
{
  const char *expl[ehisNR] = {
    "H on ND1 only",
    "H on NE2 only",
    "H on ND1 and NE2",
    "Coupled to Heme"
  };
  
  return select_res(ehisNR,resnr,hh,expl,"HISTIDINE",nrr,rr);
}

static void read_rtprename(const char *fname,FILE *fp,
                           int *nrtprename,rtprename_t **rtprename)
{
  char line[STRLEN],buf[STRLEN];
  int  n;
  rtprename_t *rr;
  int  ncol,nc;

  n  = *nrtprename;
  rr = *rtprename;

  ncol = 0;
  while(get_a_line(fp,line,STRLEN)) {
    srenew(rr,n+1);
    nc = sscanf(line,"%s %s %s %s %s %s",
                rr[n].gmx,rr[n].main,rr[n].nter,rr[n].cter,rr[n].bter,buf);
    if (ncol == 0) {
      if (nc != 2 && nc != 5) {
        gmx_fatal(FARGS,"Residue renaming database '%s' has %d columns instead of %d, %d or %d",fname,ncol,2,5);
      }
      ncol = nc;
    } else if (nc != ncol) {
        gmx_fatal(FARGS,"A line in residue renaming database '%s' has %d columns, while previous lines have %d columns",fname,nc,ncol);
    }
    
    if (nc == 2) {
      /* This file does not have special termini names, copy them from main */
      strcpy(rr[n].nter,rr[n].main);
      strcpy(rr[n].cter,rr[n].main);
      strcpy(rr[n].bter,rr[n].main);
    }

    n++;
  }

  *nrtprename = n;
  *rtprename  = rr;
}

static char *search_resrename(int nrr,rtprename_t *rr,
                              const char *name,
                              gmx_bool bStart,gmx_bool bEnd,
                              gmx_bool bCompareFFRTPname)
{
    char *nn;
    int i;

    nn = NULL;

    i = 0;
    while (i<nrr && ((!bCompareFFRTPname && strcmp(name,rr[i].gmx)  != 0) ||
                     ( bCompareFFRTPname && strcmp(name,rr[i].main) != 0)))
    {
        i++;
    }

    /* If found in the database, rename this residue's rtp building block,
     * otherwise keep the old name.
     */
    if (i < nrr)
    {
        if (bStart && bEnd)
        {
            nn = rr[i].bter;
        }
        else if (bStart)
        {
            nn = rr[i].nter;
        }
        else if (bEnd)
        {
            nn = rr[i].cter;
        }
        else
        {
            nn = rr[i].main;
        }
        if (nn[0] == '-')
        {
            gmx_fatal(FARGS,"In the chosen force field there is no residue type for '%s'%s",name,bStart ? " as a starting terminus" : (bEnd ? " as an ending terminus" : ""));
        }
    }

    return nn;
}
      

static void rename_resrtp(t_atoms *pdba,int nterpairs,int *r_start,int *r_end,
                          int nrr,rtprename_t *rr,t_symtab *symtab,
                          gmx_bool bVerbose)
{
    int  r,i,j;
    gmx_bool bStart,bEnd;
    char *nn;
    gmx_bool bFFRTPTERRNM;

    bFFRTPTERRNM = (getenv("GMX_NO_FFRTP_TER_RENAME") == NULL);

    for(r=0; r<pdba->nres; r++)
    {
        bStart = FALSE;
        bEnd   = FALSE;
        for(j=0; j<nterpairs; j++)
        {
            if (r == r_start[j])
            {
                bStart = TRUE;
            }
        }
        for(j=0; j<nterpairs; j++)
        {
            if (r == r_end[j])
            {
                bEnd = TRUE;
            }
        }

        nn = search_resrename(nrr,rr,*pdba->resinfo[r].rtp,bStart,bEnd,FALSE);

        if (bFFRTPTERRNM && nn == NULL && (bStart || bEnd))
        {
            /* This is a terminal residue, but the residue name,
             * currently stored in .rtp, is not a standard residue name,
             * but probably a force field specific rtp name.
             * Check if we need to rename it because it is terminal.
             */
            nn = search_resrename(nrr,rr,
                                  *pdba->resinfo[r].rtp,bStart,bEnd,TRUE);
        }

        if (nn != NULL && strcmp(*pdba->resinfo[r].rtp,nn) != 0)
        {
            if (bVerbose)
            {
                printf("Changing rtp entry of residue %d %s to '%s'\n",
                       pdba->resinfo[r].nr,*pdba->resinfo[r].name,nn);
            }
            pdba->resinfo[r].rtp = put_symtab(symtab,nn);
        }
    }
}

static void pdbres_to_gmxrtp(t_atoms *pdba)
{
    int i;
  
    for(i=0; (i<pdba->nres); i++)
    {
        if (pdba->resinfo[i].rtp == NULL)
        {
            pdba->resinfo[i].rtp = pdba->resinfo[i].name;
        }
    }
}

static void rename_pdbres(t_atoms *pdba,const char *oldnm,const char *newnm,
                          gmx_bool bFullCompare,t_symtab *symtab)
{
    char *resnm;
    int i;
  
    for(i=0; (i<pdba->nres); i++)
    {
        resnm = *pdba->resinfo[i].name;
        if ((bFullCompare && (gmx_strcasecmp(resnm,oldnm) == 0)) ||
            (!bFullCompare && strstr(resnm,oldnm) != NULL))
        {
            /* Rename the residue name (not the rtp name) */
            pdba->resinfo[i].name = put_symtab(symtab,newnm);
        }
    }
}

static void rename_bb(t_atoms *pdba,const char *oldnm,const char *newnm,
                      gmx_bool bFullCompare,t_symtab *symtab)
{
    char *bbnm;
    int i;
  
    for(i=0; (i<pdba->nres); i++)
    {
        /* We have not set the rtp name yes, use the residue name */
        bbnm = *pdba->resinfo[i].name;
        if ((bFullCompare && (gmx_strcasecmp(bbnm,oldnm) == 0)) ||
            (!bFullCompare && strstr(bbnm,oldnm) != NULL))
        {
            /* Change the rtp builing block name */
            pdba->resinfo[i].rtp = put_symtab(symtab,newnm);
        }
    }
}

static void rename_bbint(t_atoms *pdba,const char *oldnm,
                         const char *gettp(int,int,const rtprename_t *),
                         gmx_bool bFullCompare,
                         t_symtab *symtab,
                         int nrr,const rtprename_t *rr)
{
    int  i;
    const char *ptr;
    char *bbnm;
  
    for(i=0; i<pdba->nres; i++)
    {
        /* We have not set the rtp name yes, use the residue name */
        bbnm = *pdba->resinfo[i].name;
        if ((bFullCompare && (strcmp(bbnm,oldnm) == 0)) ||
            (!bFullCompare && strstr(bbnm,oldnm) != NULL))
        {
            ptr = gettp(i,nrr,rr);
            pdba->resinfo[i].rtp = put_symtab(symtab,ptr);
        }
    }
}

static void check_occupancy(t_atoms *atoms,const char *filename,gmx_bool bVerbose)
{
  int i,ftp;
  int nzero=0;
  int nnotone=0;
  
  ftp = fn2ftp(filename);
  if (!atoms->pdbinfo || ((ftp != efPDB) && (ftp != efBRK) && (ftp != efENT)))
    fprintf(stderr,"No occupancies in %s\n",filename);
  else {
    for(i=0; (i<atoms->nr); i++) {
      if (atoms->pdbinfo[i].occup != 1) {
        if (bVerbose)
          fprintf(stderr,"Occupancy for atom %s%d-%s is %f rather than 1\n",
                  *atoms->resinfo[atoms->atom[i].resind].name,
                  atoms->resinfo[ atoms->atom[i].resind].nr,
                  *atoms->atomname[i],
                  atoms->pdbinfo[i].occup);
        if (atoms->pdbinfo[i].occup == 0) 
          nzero++;
        else 
          nnotone++;
      }
    }
    if (nzero == atoms->nr) {
      fprintf(stderr,"All occupancy fields zero. This is probably not an X-Ray structure\n");
    } else if ((nzero > 0) || (nnotone > 0)) {
      fprintf(stderr,
              "\n"
              "WARNING: there were %d atoms with zero occupancy and %d atoms with\n"
              "         occupancy unequal to one (out of %d atoms). Check your pdb file.\n"
              "\n",
              nzero,nnotone,atoms->nr);
    } else {
      fprintf(stderr,"All occupancies are one\n");
    }
  }
}

void write_posres(char *fn,t_atoms *pdba,real fc)
{
  FILE *fp;
  int  i;
  
  fp=gmx_fio_fopen(fn,"w");
  fprintf(fp,
          "; In this topology include file, you will find position restraint\n"
          "; entries for all the heavy atoms in your original pdb file.\n"
          "; This means that all the protons which were added by pdb2gmx are\n"
          "; not restrained.\n"
          "\n"
          "[ position_restraints ]\n"
          "; %4s%6s%8s%8s%8s\n","atom","type","fx","fy","fz"
          );
  for(i=0; (i<pdba->nr); i++) {
    if (!is_hydrogen(*pdba->atomname[i]) && !is_dummymass(*pdba->atomname[i]))
      fprintf(fp,"%6d%6d  %g  %g  %g\n",i+1,1,fc,fc,fc);
  }
  gmx_fio_fclose(fp);
}

static int read_pdball(const char *inf, const char *outf,char *title,
                       t_atoms *atoms, rvec **x,
                       int *ePBC,matrix box, gmx_bool bRemoveH,
                       t_symtab *symtab,gmx_residuetype_t rt,const char *watres,
                       gmx_atomprop_t aps,gmx_bool bVerbose)
/* Read a pdb file. (containing proteins) */
{
  int  natom,new_natom,i;
  
  /* READ IT */
  printf("Reading %s...\n",inf);
  get_stx_coordnum(inf,&natom);
  init_t_atoms(atoms,natom,TRUE);
  snew(*x,natom);
  read_stx_conf(inf,title,atoms,*x,NULL,ePBC,box);
  if (fn2ftp(inf) == efPDB)
    get_pdb_atomnumber(atoms,aps);
  if (bRemoveH) {
    new_natom=0;
    for(i=0; i<atoms->nr; i++)
      if (!is_hydrogen(*atoms->atomname[i])) {
        atoms->atom[new_natom]=atoms->atom[i];
        atoms->atomname[new_natom]=atoms->atomname[i];
        atoms->pdbinfo[new_natom]=atoms->pdbinfo[i];
        copy_rvec((*x)[i],(*x)[new_natom]);
        new_natom++;
      }
    atoms->nr=new_natom;
    natom=new_natom;
  }
    
  printf("Read");
  if (title && title[0])
    printf(" '%s',",title);
  printf(" %d atoms\n",natom);
  
  /* Rename residues */
  rename_pdbres(atoms,"HOH",watres,FALSE,symtab);
  rename_pdbres(atoms,"SOL",watres,FALSE,symtab);
  rename_pdbres(atoms,"WAT",watres,FALSE,symtab);

  rename_atoms("xlateat.dat",NULL,
               atoms,symtab,NULL,TRUE,rt,TRUE,bVerbose);
  
  if (natom == 0)
    return 0;

  if (outf)
    write_sto_conf(outf,title,atoms,*x,NULL,*ePBC,box);
 
  return natom;
}

void process_chain(t_atoms *pdba, rvec *x, 
                   gmx_bool bTrpU,gmx_bool bPheU,gmx_bool bTyrU,
                   gmx_bool bLysMan,gmx_bool bAspMan,gmx_bool bGluMan,
                   gmx_bool bHisMan,gmx_bool bArgMan,gmx_bool bGlnMan,
                   real angle,real distance,t_symtab *symtab,
                   int nrr,const rtprename_t *rr)
{
  /* Rename aromatics, lys, asp and histidine */
  if (bTyrU) rename_bb(pdba,"TYR","TYRU",FALSE,symtab);
  if (bTrpU) rename_bb(pdba,"TRP","TRPU",FALSE,symtab);
  if (bPheU) rename_bb(pdba,"PHE","PHEU",FALSE,symtab);
  if (bLysMan) 
    rename_bbint(pdba,"LYS",get_lystp,FALSE,symtab,nrr,rr);
  if (bArgMan) 
    rename_bbint(pdba,"ARG",get_argtp,FALSE,symtab,nrr,rr);
  if (bGlnMan) 
    rename_bbint(pdba,"GLN",get_glntp,FALSE,symtab,nrr,rr);
  if (bAspMan) 
    rename_bbint(pdba,"ASP",get_asptp,FALSE,symtab,nrr,rr);
  else
    rename_bb(pdba,"ASPH","ASP",FALSE,symtab);
  if (bGluMan) 
    rename_bbint(pdba,"GLU",get_glutp,FALSE,symtab,nrr,rr);
  else
    rename_bb(pdba,"GLUH","GLU",FALSE,symtab);

  if (!bHisMan)
    set_histp(pdba,x,angle,distance);
  else
    rename_bbint(pdba,"HIS",get_histp,TRUE,symtab,nrr,rr);

  /* Initialize the rtp builing block names with the residue names
   * for the residues that have not been processed above.
   */
  pdbres_to_gmxrtp(pdba);

  /* Now we have all rtp names set.
   * The rtp names will conform to Gromacs naming,
   * unless the input pdb file contained one or more force field specific
   * rtp names as residue names.
   */
}

/* struct for sorting the atoms from the pdb file */
typedef struct {
  int  resnr;  /* residue number               */
  int  j;      /* database order index         */
  int  index;  /* original atom number         */
  char anm1;   /* second letter of atom name   */
  char altloc; /* alternate location indicator */
} t_pdbindex;
  
int pdbicomp(const void *a,const void *b)
{
  t_pdbindex *pa,*pb;
  int d;

  pa=(t_pdbindex *)a;
  pb=(t_pdbindex *)b;

  d = (pa->resnr - pb->resnr);
  if (d==0) {
    d = (pa->j - pb->j);
    if (d==0) {
      d = (pa->anm1 - pb->anm1);
      if (d==0)
        d = (pa->altloc - pb->altloc);
    }
  }

  return d;
}

static void sort_pdbatoms(int nrtp,t_restp restp[],t_hackblock hb[],
                          int natoms,t_atoms **pdbaptr,rvec **x,
                          t_blocka *block,char ***gnames)
{
  t_atoms *pdba,*pdbnew;
  rvec **xnew;
  int     i,j;
  t_restp *rptr;
  t_hackblock *hbr;
  t_pdbindex *pdbi;
  atom_id *a;
  char *atomnm;
  
  pdba=*pdbaptr;
  natoms=pdba->nr;
  pdbnew=NULL;
  snew(xnew,1);
  snew(pdbi, natoms);
  
  for(i=0; i<natoms; i++)
  {
      atomnm = *pdba->atomname[i];
      rptr = &restp[pdba->atom[i].resind];
      for(j=0; (j<rptr->natom); j++) 
      {
          if (gmx_strcasecmp(atomnm,*(rptr->atomname[j])) == 0) 
          {
              break;
          }
      }
      if (j==rptr->natom) 
      {
          char buf[STRLEN];
          
          sprintf(buf,
                  "Atom %s in residue %s %d was not found in rtp entry %s with %d atoms\n"
                  "while sorting atoms.\n%s",atomnm,
                  *pdba->resinfo[pdba->atom[i].resind].name,
                  pdba->resinfo[pdba->atom[i].resind].nr,
                  rptr->resname,
                  rptr->natom,
                  is_hydrogen(atomnm) ? 
                  "\nFor a hydrogen, this can be a different protonation state, or it\n"
                  "might have had a different number in the PDB file and was rebuilt\n"
                  "(it might for instance have been H3, and we only expected H1 & H2).\n"
                  "Note that hydrogens might have been added to the entry for the N-terminus.\n"
                  "Remove this hydrogen or choose a different protonation state to solve it.\n"
                  "Option -ignh will ignore all hydrogens in the input." : ".");
          gmx_fatal(FARGS,buf);
      }
      /* make shadow array to be sorted into indexgroup */
      pdbi[i].resnr  = pdba->atom[i].resind;
      pdbi[i].j      = j;
      pdbi[i].index  = i;
      pdbi[i].anm1   = atomnm[1];
      pdbi[i].altloc = pdba->pdbinfo[i].altloc;
  }
  qsort(pdbi,natoms,(size_t)sizeof(pdbi[0]),pdbicomp);
    
  /* pdba is sorted in pdbnew using the pdbi index */ 
  snew(a,natoms);
  snew(pdbnew,1);
  init_t_atoms(pdbnew,natoms,TRUE);
  snew(*xnew,natoms);
  pdbnew->nr=pdba->nr;
  pdbnew->nres=pdba->nres;
  sfree(pdbnew->resinfo);
  pdbnew->resinfo=pdba->resinfo;
  for (i=0; i<natoms; i++) {
    pdbnew->atom[i]     = pdba->atom[pdbi[i].index];
    pdbnew->atomname[i] = pdba->atomname[pdbi[i].index];
    pdbnew->pdbinfo[i]  = pdba->pdbinfo[pdbi[i].index];
    copy_rvec((*x)[pdbi[i].index],(*xnew)[i]);
     /* make indexgroup in block */
    a[i]=pdbi[i].index;
  }
  /* clean up */
  sfree(pdba->atomname);
  sfree(pdba->atom);
  sfree(pdba->pdbinfo);
  sfree(pdba);
  sfree(*x);
  /* copy the sorted pdbnew back to pdba */
  *pdbaptr=pdbnew;
  *x=*xnew;
  add_grp(block, gnames, natoms, a, "prot_sort");
  sfree(xnew);
  sfree(a);
  sfree(pdbi);
}

static int remove_duplicate_atoms(t_atoms *pdba,rvec x[],gmx_bool bVerbose)
{
  int     i,j,oldnatoms,ndel;
  t_resinfo *ri;
  
  printf("Checking for duplicate atoms....\n");
  oldnatoms    = pdba->nr;
  ndel = 0;
  /* NOTE: pdba->nr is modified inside the loop */
  for(i=1; (i < pdba->nr); i++) {
    /* compare 'i' and 'i-1', throw away 'i' if they are identical 
       this is a 'while' because multiple alternate locations can be present */
    while ( (i < pdba->nr) &&
            (pdba->atom[i-1].resind == pdba->atom[i].resind) &&
            (strcmp(*pdba->atomname[i-1],*pdba->atomname[i])==0) ) {
      ndel++;
      if (bVerbose) {
        ri = &pdba->resinfo[pdba->atom[i].resind];
        printf("deleting duplicate atom %4s  %s%4d%c",
               *pdba->atomname[i],*ri->name,ri->nr,ri->ic);
        if (ri->chainid && (ri->chainid != ' '))
          printf(" ch %c", ri->chainid);
        if (pdba->pdbinfo) {
          if (pdba->pdbinfo[i].atomnr)
            printf("  pdb nr %4d",pdba->pdbinfo[i].atomnr);
          if (pdba->pdbinfo[i].altloc && (pdba->pdbinfo[i].altloc!=' '))
            printf("  altloc %c",pdba->pdbinfo[i].altloc);
        }
        printf("\n");
      }
      pdba->nr--;
      /* We can not free, since it might be in the symtab */
      /* sfree(pdba->atomname[i]); */
      for (j=i; j < pdba->nr; j++) {
        pdba->atom[j]     = pdba->atom[j+1];
        pdba->atomname[j] = pdba->atomname[j+1];
        pdba->pdbinfo[j]  = pdba->pdbinfo[j+1];
        copy_rvec(x[j+1],x[j]);
      }
      srenew(pdba->atom,     pdba->nr);
      /* srenew(pdba->atomname, pdba->nr); */
      srenew(pdba->pdbinfo,  pdba->nr);
    }
  }
  if (pdba->nr != oldnatoms)
    printf("Now there are %d atoms. Deleted %d duplicates.\n",pdba->nr,ndel);
  
  return pdba->nr;
}

void find_nc_ter(t_atoms *pdba,int r0,int r1,int *r_start,int *r_end,gmx_residuetype_t rt)
{
    int i;
    const char *p_startrestype;
    const char *p_restype;
    int         nstartwarn,nendwarn;
    
    *r_start = -1;
    *r_end   = -1;

    nstartwarn = 0;
    nendwarn   = 0;
    
    /* Find the starting terminus (typially N or 5') */
    for(i=r0;i<r1 && *r_start==-1;i++)
    {
        gmx_residuetype_get_type(rt,*pdba->resinfo[i].name,&p_startrestype);
        if( !gmx_strcasecmp(p_startrestype,"Protein") || !gmx_strcasecmp(p_startrestype,"DNA") || !gmx_strcasecmp(p_startrestype,"RNA") )
        {
            printf("Identified residue %s%d as a starting terminus.\n",*pdba->resinfo[i].name,pdba->resinfo[i].nr);
            *r_start=i;
        }
        else 
        {            
            if(nstartwarn < 5)
            {    
                printf("Warning: Starting residue %s%d in chain not identified as Protein/RNA/DNA.\n",*pdba->resinfo[i].name,pdba->resinfo[i].nr);
            }
            if(nstartwarn == 5)
            {
                printf("More than 5 unidentified residues at start of chain - disabling further warnings.\n");
            }
            nstartwarn++;
        }
    }

    if(*r_start>=0)
    {
        /* Go through the rest of the residues, check that they are the same class, and identify the ending terminus. */
        for(i=*r_start;i<r1;i++)
        {
            gmx_residuetype_get_type(rt,*pdba->resinfo[i].name,&p_restype);
            if( !gmx_strcasecmp(p_restype,p_startrestype) && nendwarn==0)
            {
                *r_end=i;
            }
            else 
            {
                if(nendwarn < 5)
                {    
                    printf("Warning: Residue %s%d in chain has different type (%s) from starting residue %s%d (%s).\n",
                           *pdba->resinfo[i].name,pdba->resinfo[i].nr,p_restype,
                           *pdba->resinfo[*r_start].name,pdba->resinfo[*r_start].nr,p_startrestype);
                }
                if(nendwarn == 5)
                {
                    printf("More than 5 unidentified residues at end of chain - disabling further warnings.\n");
                }
                nendwarn++;                
            }
        }  
    }
    
    if(*r_end>=0)
    {
        printf("Identified residue %s%d as a ending terminus.\n",*pdba->resinfo[*r_end].name,pdba->resinfo[*r_end].nr);
    }
}


static void
modify_chain_numbers(t_atoms *       pdba,
                     const char *    chainsep)
{
    int   i;
    char  old_prev_chainid;
    char  old_this_chainid;
    int   old_prev_chainnum;
    int   old_this_chainnum;
    t_resinfo *ri;
    int   new_chainnum;
    
    enum 
    { 
        SPLIT_ID_OR_TER, 
        SPLIT_ID_AND_TER,
        SPLIT_ID_ONLY,
        SPLIT_TER_ONLY
    }
    splitting;
    
    printf("Splitting PDB chains based on ");
    splitting = SPLIT_TER_ONLY; /* keep compiler happy */
    
    /* Be a bit flexible to catch typos */
    if (!strncmp(chainsep,"id_o",4) || !strncmp(chainsep,"int",3))
    {
        /* For later interactive splitting we tentatively assign new chain numbers at either changing id or ter records */
        splitting = SPLIT_ID_OR_TER;
        printf("TER records or changing chain id.\n");
    }
    else if (!strncmp(chainsep,"id_a",4))
    {
        splitting = SPLIT_ID_AND_TER;
        printf("TER records and chain id.\n");
    }
    else if (strlen(chainsep)==2 && !strncmp(chainsep,"id",4))
    {
        splitting = SPLIT_ID_ONLY;
        printf("changing chain id only (ignoring TER records).\n");
    }
    else if (chainsep[0]=='t')
    {
        splitting = SPLIT_TER_ONLY;
        printf("TER records only (ignoring chain id).\n");
    }
    else
    {
        gmx_fatal(FARGS,"Unidentified setting for chain separation: %s\n",chainsep);
    }                                                                           
                                                                                   
    /* The default chain enumeration is based on TER records only, which is reflected in chainnum below */
    
    old_prev_chainid  = '?';
    old_prev_chainnum = -1;
    new_chainnum  = -1;
    
    for(i=0;i<pdba->nres;i++)
    {
        ri = &pdba->resinfo[i];
        old_this_chainid  = ri->chainid;
        old_this_chainnum = ri->chainnum;

        switch (splitting)
        {
            case SPLIT_ID_OR_TER:
                if(old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
                {
                    new_chainnum++;
                }
                break;
                
            case SPLIT_ID_AND_TER:
                if(old_this_chainid != old_prev_chainid && old_this_chainnum != old_prev_chainnum)
                {
                    new_chainnum++;
                }
                break;
                
            case SPLIT_ID_ONLY:
                if(old_this_chainid != old_prev_chainid)
                {
                    new_chainnum++;
                }
                break;
                
            case SPLIT_TER_ONLY:
                if(old_this_chainnum != old_prev_chainnum)
                {
                    new_chainnum++;
                }
                break;
            default:
                gmx_fatal(FARGS,"Internal inconsistency - this shouldn't happen...");
                break;
        }
        old_prev_chainid  = old_this_chainid;
        old_prev_chainnum = old_this_chainnum;
                                                                                   
        ri->chainnum = new_chainnum;        
    }
}


typedef struct {
  char chainid;
  char chainnum;
  int  start;
  int  natom;
  gmx_bool bAllWat;
  int  nterpairs;
  int  *chainstart;
} t_pdbchain;

typedef struct {
  char chainid;
  int  chainnum;
  gmx_bool bAllWat;
  int nterpairs;
  int *chainstart;
  t_hackblock **ntdb;
  t_hackblock **ctdb;
  int *r_start;
  int *r_end;
  t_atoms *pdba;
  rvec *x;
} t_chain;

int main(int argc, char *argv[])
{
  const char *desc[] = {
    "This program reads a [TT].pdb[tt] (or [TT].gro[tt]) file, reads",
    "some database files, adds hydrogens to the molecules and generates",
    "coordinates in GROMACS (GROMOS), or optionally [TT].pdb[tt], format",
    "and a topology in GROMACS format.",
    "These files can subsequently be processed to generate a run input file.",
    "[PAR]",
    "[TT]pdb2gmx[tt] will search for force fields by looking for",
    "a [TT]forcefield.itp[tt] file in subdirectories [TT]<forcefield>.ff[tt]",
    "of the current working directory and of the GROMACS library directory",
    "as inferred from the path of the binary or the [TT]GMXLIB[tt] environment",
    "variable.",
    "By default the forcefield selection is interactive,",
    "but you can use the [TT]-ff[tt] option to specify one of the short names",
    "in the list on the command line instead. In that case [TT]pdb2gmx[tt] just looks",
    "for the corresponding [TT]<forcefield>.ff[tt] directory.",
    "[PAR]",
    "After choosing a force field, all files will be read only from",
    "the corresponding force field directory.",
    "If you want to modify or add a residue types, you can copy the force",
    "field directory from the GROMACS library directory to your current",
    "working directory. If you want to add new protein residue types,",
    "you will need to modify [TT]residuetypes.dat[tt] in the library directory",
    "or copy the whole library directory to a local directory and set",
    "the environment variable [TT]GMXLIB[tt] to the name of that directory.",
    "Check Chapter 5 of the manual for more information about file formats.",
    "[PAR]",
    
    "Note that a [TT].pdb[tt] file is nothing more than a file format, and it",
    "need not necessarily contain a protein structure. Every kind of",
    "molecule for which there is support in the database can be converted.",
    "If there is no support in the database, you can add it yourself.[PAR]",
    
    "The program has limited intelligence, it reads a number of database",
    "files, that allow it to make special bonds (Cys-Cys, Heme-His, etc.),",
    "if necessary this can be done manually. The program can prompt the",
    "user to select which kind of LYS, ASP, GLU, CYS or HIS residue is",
    "desired. For Lys the choice is between neutral (two protons on NZ) or",
    "protonated (three protons, default), for Asp and Glu unprotonated",
    "(default) or protonated, for His the proton can be either on ND1,",
    "on NE2 or on both. By default these selections are done automatically.",
    "For His, this is based on an optimal hydrogen bonding",
    "conformation. Hydrogen bonds are defined based on a simple geometric",
    "criterion, specified by the maximum hydrogen-donor-acceptor angle",
    "and donor-acceptor distance, which are set by [TT]-angle[tt] and",
    "[TT]-dist[tt] respectively.[PAR]",
     
    "The protonation state of N- and C-termini can be chosen interactively",
    "with the [TT]-ter[tt] flag.  Default termini are ionized (NH3+ and COO-),",
    "respectively.  Some force fields support zwitterionic forms of isolated",
    "amino acids, but for polypeptides these options should NOT be selected.[PAR]",
 
    "The separation of chains is not entirely trivial since the markup",
    "in user-generated PDB files frequently varies and sometimes it",
    "is desirable to merge entries across a TER record, for instance",
    "if you want a disulfide bridge or distance restraints between",
    "two protein chains or if you have a HEME group bound to a protein.",
    "In such cases multiple chains should be contained in a single",
    "[TT]moleculetype[tt] definition.",
    "To handle this, [TT]pdb2gmx[tt] has an option [TT]-chainsep[tt] so you can",
    "choose whether a new chain should start when we find a TER record,",
    "when the chain id changes, combinations of either or both of these",
    "or fully interactively.[PAR]",
    
    "[TT]pdb2gmx[tt] will also check the occupancy field of the [TT].pdb[tt] file.",
    "If any of the occupancies are not one, indicating that the atom is",
    "not resolved well in the structure, a warning message is issued.",
    "When a [TT].pdb[tt] file does not originate from an X-ray structure determination",
    "all occupancy fields may be zero. Either way, it is up to the user",
    "to verify the correctness of the input data (read the article!).[PAR]", 
    
    "During processing the atoms will be reordered according to GROMACS",
    "conventions. With [TT]-n[tt] an index file can be generated that",
    "contains one group reordered in the same way. This allows you to",
    "convert a GROMOS trajectory and coordinate file to GROMOS. There is",
    "one limitation: reordering is done after the hydrogens are stripped",
    "from the input and before new hydrogens are added. This means that",
    "you should not use [TT]-ignh[tt].[PAR]",

    "The [TT].gro[tt] and [TT].g96[tt] file formats do not support chain",
    "identifiers. Therefore it is useful to enter a [TT].pdb[tt] file name at",
    "the [TT]-o[tt] option when you want to convert a multi-chain [TT].pdb[tt] file.",
    "[PAR]",
    
    "The option [TT]-vsite[tt] removes hydrogen and fast improper dihedral",
    "motions. Angular and out-of-plane motions can be removed by changing",
    "hydrogens into virtual sites and fixing angles, which fixes their",
    "position relative to neighboring atoms. Additionally, all atoms in the",
    "aromatic rings of the standard amino acids (i.e. PHE, TRP, TYR and HIS)",
    "can be converted into virtual sites, eliminating the fast improper dihedral",
    "fluctuations in these rings. [BB]Note[bb] that in this case all other hydrogen",
    "atoms are also converted to virtual sites. The mass of all atoms that are",
    "converted into virtual sites, is added to the heavy atoms.[PAR]",
    "Also slowing down of dihedral motion can be done with [TT]-heavyh[tt]",
    "done by increasing the hydrogen-mass by a factor of 4. This is also",
    "done for water hydrogens to slow down the rotational motion of water.",
    "The increase in mass of the hydrogens is subtracted from the bonded",
    "(heavy) atom so that the total mass of the system remains the same."
  };

  
  FILE       *fp,*top_file,*top_file2,*itp_file=NULL;
  int        natom,nres;
  t_atoms    pdba_all,*pdba;
  t_atoms    *atoms;
  t_resinfo  *ri;
  t_blocka   *block;
  int        chain,nch,maxch,nwaterchain;
  t_pdbchain *pdb_ch;
  t_chain    *chains,*cc;
  char       select[STRLEN];
  int        nincl,nmol;
  char       **incls;
  t_mols     *mols;
  char       **gnames;
  int        ePBC;
  matrix     box;
  rvec       box_space;
  int        i,j,k,l,nrtp;
  int        *swap_index,si;
  t_restp    *restp;
  t_hackblock *ah;
  t_symtab   symtab;
  gpp_atomtype_t atype;
  gmx_residuetype_t rt;
  const char *top_fn;
  char       fn[256],itp_fn[STRLEN],posre_fn[STRLEN],buf_fn[STRLEN];
  char       molname[STRLEN],title[STRLEN],quote[STRLEN],generator[STRLEN];
  char       *c,forcefield[STRLEN],ffdir[STRLEN];
  char       ffname[STRLEN],suffix[STRLEN],buf[STRLEN];
  char       *watermodel;
  const char *watres;
  int        nrtpf;
  char       **rtpf;
  char       rtp[STRLEN];
  int        nrrn;
  char       **rrn;
  int        nrtprename,naa;
  rtprename_t *rtprename=NULL;
  int        nah,nNtdb,nCtdb,ntdblist;
  t_hackblock *ntdb,*ctdb,**tdblist;
  int        nssbonds;
  t_ssbond   *ssbonds;
  rvec       *pdbx,*x;
  gmx_bool       bVsites=FALSE,bWat,bPrevWat=FALSE,bITP,bVsiteAromatics=FALSE,bMerge;
  real       mHmult=0;
  t_hackblock *hb_chain;
  t_restp    *restp_chain;
  output_env_t oenv;
  const char *p_restype;
  int        rc;
  int           this_atomnum;
  int           prev_atomnum;
  const char *  prev_atomname;
  const char *  this_atomname;
  const char *  prev_resname;
  const char *  this_resname;
  int           prev_resnum;
  int           this_resnum;
  char          prev_chainid;
  char          this_chainid;
  int           prev_chainnumber;
  int           this_chainnumber;
  int           nid_used;
  int           this_chainstart;
  int           prev_chainstart;
  gmx_bool      bMerged;
    
  gmx_atomprop_t aps;
  
  t_filenm   fnm[] = { 
    { efSTX, "-f", "eiwit.pdb", ffREAD  },
    { efSTO, "-o", "conf",      ffWRITE },
    { efTOP, NULL, NULL,        ffWRITE },
    { efITP, "-i", "posre",     ffWRITE },
    { efNDX, "-n", "clean",     ffOPTWR },
    { efSTO, "-q", "clean.pdb", ffOPTWR }
  };
#define NFILE asize(fnm)
 

  /* Command line arguments must be static */
  static gmx_bool bNewRTP=FALSE;
  static gmx_bool bInter=FALSE, bCysMan=FALSE; 
  static gmx_bool bLysMan=FALSE, bAspMan=FALSE, bGluMan=FALSE, bHisMan=FALSE;
  static gmx_bool bGlnMan=FALSE, bArgMan=FALSE;
  static gmx_bool bTerMan=FALSE, bUnA=FALSE, bHeavyH;
  static gmx_bool bSort=TRUE, bAllowMissing=FALSE, bRemoveH=FALSE;
  static gmx_bool bDeuterate=FALSE,bVerbose=FALSE,bChargeGroups=TRUE,bCmap=TRUE;
  static gmx_bool bRenumRes=FALSE,bRTPresname=FALSE;
  static real angle=135.0, distance=0.3,posre_fc=1000;
  static real long_bond_dist=0.25, short_bond_dist=0.05;
  static const char *vsitestr[] = { NULL, "none", "hydrogens", "aromatics", NULL };
  static const char *watstr[] = { NULL, "select", "none", "spc", "spce", "tip3p", "tip4p", "tip5p", NULL };
  static const char *chainsep[] = { NULL, "id_or_ter", "id_and_ter", "ter", "id", "interactive", NULL };
  static const char *ff = "select";

  t_pargs pa[] = {
    { "-newrtp", FALSE, etBOOL, {&bNewRTP},
      "HIDDENWrite the residue database in new format to [TT]new.rtp[tt]"},
    { "-lb",     FALSE, etREAL, {&long_bond_dist},
      "HIDDENLong bond warning distance" },
    { "-sb",     FALSE, etREAL, {&short_bond_dist},
      "HIDDENShort bond warning distance" },
    { "-chainsep", FALSE, etENUM, {chainsep},
      "Condition in PDB files when a new chain and molecule_type should be started" },
    { "-ff",     FALSE, etSTR,  {&ff},
      "Force field, interactive by default. Use [TT]-h[tt] for information." },
    { "-water",  FALSE, etENUM, {watstr},
      "Water model to use" },
    { "-inter",  FALSE, etBOOL, {&bInter},
      "Set the next 8 options to interactive"},
    { "-ss",     FALSE, etBOOL, {&bCysMan}, 
      "Interactive SS bridge selection" },
    { "-ter",    FALSE, etBOOL, {&bTerMan}, 
      "Interactive termini selection, iso charged" },
    { "-lys",    FALSE, etBOOL, {&bLysMan}, 
      "Interactive lysine selection, iso charged" },
    { "-arg",    FALSE, etBOOL, {&bArgMan}, 
      "Interactive arginine selection, iso charged" },
    { "-asp",    FALSE, etBOOL, {&bAspMan}, 
      "Interactive aspartic Acid selection, iso charged" },
    { "-glu",    FALSE, etBOOL, {&bGluMan}, 
      "Interactive glutamic Acid selection, iso charged" },
    { "-gln",    FALSE, etBOOL, {&bGlnMan}, 
      "Interactive glutamine selection, iso neutral" },
    { "-his",    FALSE, etBOOL, {&bHisMan},
      "Interactive histidine selection, iso checking H-bonds" },
    { "-angle",  FALSE, etREAL, {&angle}, 
      "Minimum hydrogen-donor-acceptor angle for a H-bond (degrees)" },
    { "-dist",   FALSE, etREAL, {&distance},
      "Maximum donor-acceptor distance for a H-bond (nm)" },
    { "-una",    FALSE, etBOOL, {&bUnA}, 
      "Select aromatic rings with united CH atoms on phenylalanine, "
      "tryptophane and tyrosine" },
    { "-sort",   FALSE, etBOOL, {&bSort}, 
      "HIDDENSort the residues according to database, turning this off is dangerous as charge groups might be broken in parts" },
    { "-ignh",   FALSE, etBOOL, {&bRemoveH}, 
      "Ignore hydrogen atoms that are in the coordinate file" },
    { "-missing",FALSE, etBOOL, {&bAllowMissing}, 
      "Continue when atoms are missing, dangerous" },
    { "-v",      FALSE, etBOOL, {&bVerbose}, 
      "Be slightly more verbose in messages" },
    { "-posrefc",FALSE, etREAL, {&posre_fc},
      "Force constant for position restraints" },
    { "-vsite",  FALSE, etENUM, {vsitestr}, 
      "Convert atoms to virtual sites" },
    { "-heavyh", FALSE, etBOOL, {&bHeavyH},
      "Make hydrogen atoms heavy" },
    { "-deuterate", FALSE, etBOOL, {&bDeuterate},
      "Change the mass of hydrogens to 2 amu" },
    { "-chargegrp", TRUE, etBOOL, {&bChargeGroups},
      "Use charge groups in the [TT].rtp[tt] file"  },
    { "-cmap", TRUE, etBOOL, {&bCmap},
      "Use cmap torsions (if enabled in the [TT].rtp[tt] file)"  },
    { "-renum", TRUE, etBOOL, {&bRenumRes},
      "Renumber the residues consecutively in the output"  },
    { "-rtpres", TRUE, etBOOL, {&bRTPresname},
      "Use [TT].rtp[tt] entry names as residue names"  }
  };
#define NPARGS asize(pa)
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,asize(desc),desc,
                    0,NULL,&oenv);

  /* Force field selection, interactive or direct */
  choose_ff(strcmp(ff,"select") == 0 ? NULL : ff,
            forcefield,sizeof(forcefield),
            ffdir,sizeof(ffdir));

  if (strlen(forcefield) > 0) {
    strcpy(ffname,forcefield);
    ffname[0] = toupper(ffname[0]);
  } else {
    gmx_fatal(FARGS,"Empty forcefield string");
  }
  
  printf("\nUsing the %s force field in directory %s\n\n",
         ffname,ffdir);
    
  choose_watermodel(watstr[0],ffdir,&watermodel);

  if (bInter) {
    /* if anything changes here, also change description of -inter */
    bCysMan = TRUE;
    bTerMan = TRUE;
    bLysMan = TRUE;
    bArgMan = TRUE;
    bAspMan = TRUE;
    bGluMan = TRUE;
    bGlnMan = TRUE;
    bHisMan = TRUE;
  }
  
  if (bHeavyH)
    mHmult=4.0;
  else if (bDeuterate)
    mHmult=2.0;
  else
    mHmult=1.0;
  
  switch(vsitestr[0][0]) {
  case 'n': /* none */
    bVsites=FALSE;
    bVsiteAromatics=FALSE;
    break;
  case 'h': /* hydrogens */
    bVsites=TRUE;
    bVsiteAromatics=FALSE;
    break;
  case 'a': /* aromatics */
    bVsites=TRUE;
    bVsiteAromatics=TRUE;
    break;
  default:
    gmx_fatal(FARGS,"DEATH HORROR in $s (%d): vsitestr[0]='%s'",
                __FILE__,__LINE__,vsitestr[0]);
  }/* end switch */
  
  /* Open the symbol table */
  open_symtab(&symtab);

  /* Residue type database */  
  gmx_residuetype_init(&rt);
  
  /* Read residue renaming database(s), if present */
  nrrn = fflib_search_file_end(ffdir,".r2b",FALSE,&rrn);
    
  nrtprename = 0;
  rtprename  = NULL;
  for(i=0; i<nrrn; i++) {
    fp = fflib_open(rrn[i]);
    read_rtprename(rrn[i],fp,&nrtprename,&rtprename);
    fclose(fp);
    sfree(rrn[i]);
  }
  sfree(rrn);

  /* Add all alternative names from the residue renaming database to the list of recognized amino/nucleic acids. */
  naa=0;
  for(i=0;i<nrtprename;i++)
  {
      rc=gmx_residuetype_get_type(rt,rtprename[i].gmx,&p_restype);

      /* Only add names if the 'standard' gromacs/iupac base name was found */
      if(rc==0)
      {
          gmx_residuetype_add(rt,rtprename[i].main,p_restype);
          gmx_residuetype_add(rt,rtprename[i].nter,p_restype);
          gmx_residuetype_add(rt,rtprename[i].cter,p_restype);
          gmx_residuetype_add(rt,rtprename[i].bter,p_restype);
      }          
  }
    
  clear_mat(box);
  if (watermodel != NULL && (strstr(watermodel,"4p") ||
                             strstr(watermodel,"4P"))) {
    watres = "HO4";
  } else if (watermodel != NULL && (strstr(watermodel,"5p") ||
                                    strstr(watermodel,"5P"))) {
    watres = "HO5";
  } else {
    watres = "HOH";
  }
    
  aps = gmx_atomprop_init();
  natom = read_pdball(opt2fn("-f",NFILE,fnm),opt2fn_null("-q",NFILE,fnm),title,
                      &pdba_all,&pdbx,&ePBC,box,bRemoveH,&symtab,rt,watres,
                      aps,bVerbose);
  
  if (natom==0)
    gmx_fatal(FARGS,"No atoms found in pdb file %s\n",opt2fn("-f",NFILE,fnm));

  printf("Analyzing pdb file\n");
  nch=0;
  maxch=0;
  nwaterchain=0;
    
  modify_chain_numbers(&pdba_all,chainsep[0]);

    
  bMerge = !strncmp(chainsep[0],"int",3);
    
  this_atomname       = NULL;
  this_atomnum        = -1;
  this_resname        = NULL;
  this_resnum         = -1;
  this_chainid        = '?';
  this_chainnumber    = -1;
  this_chainstart     = 0;
  /* Keep the compiler happy */
  prev_chainstart     = 0;
    
  pdb_ch=NULL;
  bMerged = FALSE;
  for (i=0; (i<natom); i++) 
  {
      ri = &pdba_all.resinfo[pdba_all.atom[i].resind];

      prev_atomname      = this_atomname;
      prev_atomnum       = this_atomnum;
      prev_resname       = this_resname;
      prev_resnum        = this_resnum;
      prev_chainid       = this_chainid;
      prev_chainnumber   = this_chainnumber;
      if (!bMerged)
      {
          prev_chainstart    = this_chainstart;
      }
      
      this_atomname      = *pdba_all.atomname[i];
      this_atomnum       = (pdba_all.pdbinfo != NULL) ? pdba_all.pdbinfo[i].atomnr : i+1;
      this_resname       = *ri->name;
      this_resnum        = ri->nr;
      this_chainid       = ri->chainid;
      this_chainnumber   = ri->chainnum;
      
      bWat = gmx_strcasecmp(*ri->name,watres) == 0;
      if ((i == 0) || (this_chainnumber != prev_chainnumber) || (bWat != bPrevWat)) 
      {
          this_chainstart = pdba_all.atom[i].resind;
          if (bMerge && i>0 && !bWat) 
          {
              printf("Merge chain ending with residue %s%d (chain id '%c', atom %d %s) with\n"
                     "chain starting with residue %s%d (chain id '%c', atom %d %s)? [n/y]\n",
                     prev_resname,prev_resnum,prev_chainid,prev_atomnum,prev_atomname,
                     this_resname,this_resnum,this_chainid,this_atomnum,this_atomname);
              
              if(NULL==fgets(select,STRLEN-1,stdin))
              {
                  gmx_fatal(FARGS,"Error reading from stdin");
              }
          } 
          else
          {
              select[0] = 'n';
          }
          
          bMerged = (select[0] == 'y');
          if (bMerged) 
          {
              pdb_ch[nch-1].chainstart[pdb_ch[nch-1].nterpairs] = 
                  pdba_all.atom[i].resind - prev_chainstart;
              pdb_ch[nch-1].nterpairs++;
              srenew(pdb_ch[nch-1].chainstart,pdb_ch[nch-1].nterpairs+1);
          }
          else 
          {
              /* set natom for previous chain */
              if (nch > 0)
              {
                  pdb_ch[nch-1].natom=i-pdb_ch[nch-1].start;
              }
              if (bWat)
              {
                  nwaterchain++;
                  ri->chainid = ' ';
              }
              /* check if chain identifier was used before */
              for (j=0; (j<nch); j++) 
              {
                  if (pdb_ch[j].chainid != ' ' && pdb_ch[j].chainid == ri->chainid) 
                  {
                      printf("WARNING: Chain identifier '%c' is used in two non-sequential blocks.\n"
                             "They will be treated as separate chains unless you reorder your file.\n",
                             ri->chainid);
                  }
              }
              if (nch == maxch)
              {
                  maxch += 16;
                  srenew(pdb_ch,maxch);
              }
              pdb_ch[nch].chainid = ri->chainid;
              pdb_ch[nch].chainnum = ri->chainnum; 
              pdb_ch[nch].start=i;
              pdb_ch[nch].bAllWat=bWat;
              if (bWat)
                  pdb_ch[nch].nterpairs=0;
              else
                  pdb_ch[nch].nterpairs=1;
              snew(pdb_ch[nch].chainstart,pdb_ch[nch].nterpairs+1);
              /* modified [nch] to [0] below */
              pdb_ch[nch].chainstart[0]=0;
              nch++;
          }
      }
      bPrevWat=bWat;
  }
  pdb_ch[nch-1].natom=natom-pdb_ch[nch-1].start;
  
  /* set all the water blocks at the end of the chain */
  snew(swap_index,nch);
  j=0;
  for(i=0; i<nch; i++)
    if (!pdb_ch[i].bAllWat) {
      swap_index[j]=i;
      j++;
    }
  for(i=0; i<nch; i++)
    if (pdb_ch[i].bAllWat) {
      swap_index[j]=i;
      j++;
    }
  if (nwaterchain>1)
    printf("Moved all the water blocks to the end\n");

  snew(chains,nch);
  /* copy pdb data and x for all chains */
  for (i=0; (i<nch); i++) {
    si=swap_index[i];
    chains[i].chainid = pdb_ch[si].chainid;
    chains[i].chainnum = pdb_ch[si].chainnum;
    chains[i].bAllWat = pdb_ch[si].bAllWat;
    chains[i].nterpairs = pdb_ch[si].nterpairs;
    chains[i].chainstart = pdb_ch[si].chainstart;
    snew(chains[i].ntdb,pdb_ch[si].nterpairs);
    snew(chains[i].ctdb,pdb_ch[si].nterpairs);
    snew(chains[i].r_start,pdb_ch[si].nterpairs);
    snew(chains[i].r_end,pdb_ch[si].nterpairs);
      
    snew(chains[i].pdba,1);
    init_t_atoms(chains[i].pdba,pdb_ch[si].natom,TRUE);
    snew(chains[i].x,chains[i].pdba->nr);
    for (j=0; j<chains[i].pdba->nr; j++) {
      chains[i].pdba->atom[j] = pdba_all.atom[pdb_ch[si].start+j];
      snew(chains[i].pdba->atomname[j],1);
      *chains[i].pdba->atomname[j] = 
        strdup(*pdba_all.atomname[pdb_ch[si].start+j]);
      chains[i].pdba->pdbinfo[j] = pdba_all.pdbinfo[pdb_ch[si].start+j];
      copy_rvec(pdbx[pdb_ch[si].start+j],chains[i].x[j]);
    }
    /* Re-index the residues assuming that the indices are continuous */
    k    = chains[i].pdba->atom[0].resind;
    nres = chains[i].pdba->atom[chains[i].pdba->nr-1].resind - k + 1;
    chains[i].pdba->nres = nres;
    for(j=0; j < chains[i].pdba->nr; j++) {
      chains[i].pdba->atom[j].resind -= k;
    }
    srenew(chains[i].pdba->resinfo,nres);
    for(j=0; j<nres; j++) {
      chains[i].pdba->resinfo[j] = pdba_all.resinfo[k+j];
      snew(chains[i].pdba->resinfo[j].name,1);
      *chains[i].pdba->resinfo[j].name = strdup(*pdba_all.resinfo[k+j].name);
      /* make all chain identifiers equal to that of the chain */
      chains[i].pdba->resinfo[j].chainid = pdb_ch[si].chainid;
    }
  }

  if (bMerge)
    printf("\nMerged %d chains into one molecule definition\n\n",
           pdb_ch[0].nterpairs);

  printf("There are %d chains and %d blocks of water and "
         "%d residues with %d atoms\n",
         nch-nwaterchain,nwaterchain,
         pdba_all.resinfo[pdba_all.atom[natom-1].resind].nr,natom);
          
  printf("\n  %5s  %4s %6s\n","chain","#res","#atoms");
  for (i=0; (i<nch); i++)
    printf("  %d '%c' %5d %6d  %s\n",
           i+1, chains[i].chainid ? chains[i].chainid:'-',
           chains[i].pdba->nres, chains[i].pdba->nr,
           chains[i].bAllWat ? "(only water)":"");
  printf("\n");
  
  check_occupancy(&pdba_all,opt2fn("-f",NFILE,fnm),bVerbose);
  
  /* Read atomtypes... */
  atype = read_atype(ffdir,&symtab);
  
  /* read residue database */
  printf("Reading residue database... (%s)\n",forcefield);
  nrtpf = fflib_search_file_end(ffdir,".rtp",TRUE,&rtpf);
  nrtp  = 0;
  restp = NULL;
  for(i=0; i<nrtpf; i++) {
    read_resall(rtpf[i],&nrtp,&restp,atype,&symtab,FALSE);
    sfree(rtpf[i]);
  }
  sfree(rtpf);
  if (bNewRTP) {
    /* Not correct with multiple rtp input files with different bonded types */
    fp=gmx_fio_fopen("new.rtp","w");
    print_resall(fp,nrtp,restp,atype);
    gmx_fio_fclose(fp);
  }
    
  /* read hydrogen database */
  nah = read_h_db(ffdir,&ah);
  
  /* Read Termini database... */
  nNtdb=read_ter_db(ffdir,'n',&ntdb,atype);
  nCtdb=read_ter_db(ffdir,'c',&ctdb,atype);
  
  top_fn=ftp2fn(efTOP,NFILE,fnm);
  top_file=gmx_fio_fopen(top_fn,"w");

  sprintf(generator,"%s - %s",ShortProgram(), GromacsVersion() );

  print_top_header(top_file,top_fn,generator,FALSE,ffdir,mHmult);

  nincl=0;
  nmol=0;
  incls=NULL;
  mols=NULL;
  nres=0;
  for(chain=0; (chain<nch); chain++) {
    cc = &(chains[chain]);

    /* set pdba, natom and nres to the current chain */
    pdba =cc->pdba;
    x    =cc->x;
    natom=cc->pdba->nr;
    nres =cc->pdba->nres;
    
    if (cc->chainid && ( cc->chainid != ' ' ) )
      printf("Processing chain %d '%c' (%d atoms, %d residues)\n",
              chain+1,cc->chainid,natom,nres);
    else
      printf("Processing chain %d (%d atoms, %d residues)\n",
              chain+1,natom,nres);
      
    process_chain(pdba,x,bUnA,bUnA,bUnA,bLysMan,bAspMan,bGluMan,
                  bHisMan,bArgMan,bGlnMan,angle,distance,&symtab,
                  nrtprename,rtprename);
      
        cc->chainstart[cc->nterpairs] = pdba->nres;
        j = 0;
        for(i=0; i<cc->nterpairs; i++)
        {
            find_nc_ter(pdba,cc->chainstart[i],cc->chainstart[i+1],
                        &(cc->r_start[j]),&(cc->r_end[j]),rt);    
      
            if (cc->r_start[j] >= 0 && cc->r_end[j] >= 0)
            {
                j++;
            }
        }
        cc->nterpairs = j;
        if (cc->nterpairs == 0)
        {
            printf("Problem with chain definition, or missing terminal residues.\n"
                   "This chain does not appear to contain a recognized chain molecule.\n"
                   "If this is incorrect, you can edit residuetypes.dat to modify the behavior.\n");
        }

    /* Check for disulfides and other special bonds */
    nssbonds = mk_specbonds(pdba,x,bCysMan,&ssbonds,bVerbose);

    if (nrtprename > 0) {        
      rename_resrtp(pdba,cc->nterpairs,cc->r_start,cc->r_end,nrtprename,rtprename,
                    &symtab,bVerbose);
    }
    
    if (debug) {
      if (nch==1) {
        sprintf(fn,"chain.pdb");
      } else {
        sprintf(fn,"chain_%c%d.pdb",cc->chainid,cc->chainnum);
      }
      write_sto_conf(fn,title,pdba,x,NULL,ePBC,box);
    }

      
    for(i=0; i<cc->nterpairs; i++) 
    {
        
        /* Set termini.
         * We first apply a filter so we only have the
         * termini that can be applied to the residue in question
         * (or a generic terminus if no-residue specific is available).
         */
        /* First the N terminus */
        if (nNtdb > 0) 
        {
            tdblist = filter_ter(nrtp,restp,nNtdb,ntdb,
                                 *pdba->resinfo[cc->r_start[i]].name,
                                 *pdba->resinfo[cc->r_start[i]].rtp,
                                 &ntdblist);
            if(ntdblist==0)
            {
                printf("No suitable end (N or 5') terminus found in database - assuming this residue\n"
                       "is already in a terminus-specific form and skipping terminus selection.\n");
                cc->ntdb[i]=NULL;
            }
            else 
            {
                if(bTerMan && ntdblist>1)
                {
                    sprintf(select,"Select start terminus type for %s-%d",
                            *pdba->resinfo[cc->r_start[i]].name,
                            pdba->resinfo[cc->r_start[i]].nr);
                    cc->ntdb[i] = choose_ter(ntdblist,tdblist,select);
                }
                else
                {
                    cc->ntdb[i] = tdblist[0];
                }
                
                printf("Start terminus %s-%d: %s\n",
                       *pdba->resinfo[cc->r_start[i]].name,
                       pdba->resinfo[cc->r_start[i]].nr,
                       (cc->ntdb[i])->name);
                sfree(tdblist);
            }
        }
        else 
        {
            cc->ntdb[i] = NULL;
        }
        
        /* And the C terminus */
        if (nCtdb > 0)
        {
            tdblist = filter_ter(nrtp,restp,nCtdb,ctdb,
                                 *pdba->resinfo[cc->r_end[i]].name,
                                 *pdba->resinfo[cc->r_end[i]].rtp,
                                 &ntdblist);
            if(ntdblist==0)
            {
                printf("No suitable end (C or 3') terminus found in database - assuming this residue\n"
                       "is already in a terminus-specific form and skipping terminus selection.\n");
                cc->ctdb[i] = NULL;
            }
            else 
            {
                if(bTerMan && ntdblist>1)
                {
                    sprintf(select,"Select end terminus type for %s-%d",
                            *pdba->resinfo[cc->r_end[i]].name,
                            pdba->resinfo[cc->r_end[i]].nr);
                    cc->ctdb[i] = choose_ter(ntdblist,tdblist,select);
                }
                else
                {
                    cc->ctdb[i] = tdblist[0];
                }
                printf("End terminus %s-%d: %s\n",
                       *pdba->resinfo[cc->r_end[i]].name,
                       pdba->resinfo[cc->r_end[i]].nr,
                       (cc->ctdb[i])->name);
                sfree(tdblist);
            }
        }
        else 
        {
            cc->ctdb[i] = NULL;
        }
    }
    /* lookup hackblocks and rtp for all residues */
    get_hackblocks_rtp(&hb_chain, &restp_chain,
                       nrtp, restp, pdba->nres, pdba->resinfo, 
                       cc->nterpairs, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end);
    /* ideally, now we would not need the rtp itself anymore, but do 
     everything using the hb and restp arrays. Unfortunately, that 
     requires some re-thinking of code in gen_vsite.c, which I won't 
     do now :( AF 26-7-99 */

    rename_atoms(NULL,ffdir,
                 pdba,&symtab,restp_chain,FALSE,rt,FALSE,bVerbose);

    match_atomnames_with_rtp(restp_chain,hb_chain,pdba,x,bVerbose);

    if (bSort) {
      block = new_blocka();
      snew(gnames,1);
      sort_pdbatoms(pdba->nres,restp_chain,hb_chain,
                    natom,&pdba,&x,block,&gnames);
      natom = remove_duplicate_atoms(pdba,x,bVerbose);
      if (ftp2bSet(efNDX,NFILE,fnm)) {
        if (bRemoveH) {
          fprintf(stderr,"WARNING: with the -remh option the generated "
                  "index file (%s) might be useless\n"
                  "(the index file is generated before hydrogens are added)",
                  ftp2fn(efNDX,NFILE,fnm));
        }
        write_index(ftp2fn(efNDX,NFILE,fnm),block,gnames);
      }
      for(i=0; i < block->nr; i++)
        sfree(gnames[i]);
      sfree(gnames);
      done_blocka(block);
    } else {
      fprintf(stderr,"WARNING: "
              "without sorting no check for duplicate atoms can be done\n");
    }

    /* Generate Hydrogen atoms (and termini) in the sequence */
    natom=add_h(&pdba,&x,nah,ah,
                cc->nterpairs,cc->ntdb,cc->ctdb,cc->r_start,cc->r_end,bAllowMissing,
                NULL,NULL,TRUE,FALSE);
    printf("Now there are %d residues with %d atoms\n",
           pdba->nres,pdba->nr);
    if (debug) write_pdbfile(debug,title,pdba,x,ePBC,box,' ',0,NULL,TRUE);

    if (debug)
      for(i=0; (i<natom); i++)
        fprintf(debug,"Res %s%d atom %d %s\n",
                *(pdba->resinfo[pdba->atom[i].resind].name),
                pdba->resinfo[pdba->atom[i].resind].nr,i+1,*pdba->atomname[i]);
    
    strcpy(posre_fn,ftp2fn(efITP,NFILE,fnm));
    
    /* make up molecule name(s) */

      k = (cc->nterpairs>0 && cc->r_start[0]>=0) ? cc->r_start[0] : 0;
            
    gmx_residuetype_get_type(rt,*pdba->resinfo[k].name,&p_restype);
      
    suffix[0]='\0';
      
    if (cc->bAllWat) 
    {
        sprintf(molname,"Water");
    } 
    else
    {
        this_chainid = cc->chainid;
        
        /* Add the chain id if we have one */
        if(this_chainid != ' ')
        {
            sprintf(buf,"_chain_%c",this_chainid);
            strcat(suffix,buf);
        }

        /* Check if there have been previous chains with the same id */
        nid_used = 0;
        for(k=0;k<chain;k++)
        {
            if(cc->chainid == chains[k].chainid)
            {
                nid_used++;
            }
        }
        /* Add the number for this chain identifier if there are multiple copies */
        if(nid_used>0)
        {
            
            sprintf(buf,"%d",nid_used+1);
            strcat(suffix,buf);
        }

        if(strlen(suffix)>0)
        {
            sprintf(molname,"%s%s",p_restype,suffix);
        }
        else
        {
            strcpy(molname,p_restype);
        }
    }
      
    if ((nch-nwaterchain>1) && !cc->bAllWat) {
      bITP=TRUE;
      strcpy(itp_fn,top_fn);
      printf("Chain time...\n");
      c=strrchr(itp_fn,'.');
      sprintf(c,"_%s.itp",molname);
      c=strrchr(posre_fn,'.');
      sprintf(c,"_%s.itp",molname);
      if (strcmp(itp_fn,posre_fn) == 0) {
        strcpy(buf_fn,posre_fn);
        c  = strrchr(buf_fn,'.');
        *c = '\0';
        sprintf(posre_fn,"%s_pr.itp",buf_fn);
      }
      
      nincl++;
      srenew(incls,nincl);
      incls[nincl-1]=strdup(itp_fn);
      itp_file=gmx_fio_fopen(itp_fn,"w");
    } else
      bITP=FALSE;

    srenew(mols,nmol+1);
    if (cc->bAllWat) {
      mols[nmol].name = strdup("SOL");
      mols[nmol].nr   = pdba->nres;
    } else {
      mols[nmol].name = strdup(molname);
      mols[nmol].nr   = 1;
    }
    nmol++;

    if (bITP)
      print_top_comment(itp_file,itp_fn,generator,ffdir,TRUE);

    if (cc->bAllWat)
      top_file2=NULL;
    else
      if (bITP)
        top_file2=itp_file;
      else
        top_file2=top_file;

    pdb2top(top_file2,posre_fn,molname,pdba,&x,atype,&symtab,
            nrtp,restp,
            restp_chain,hb_chain,
            cc->nterpairs,cc->ntdb,cc->ctdb,bAllowMissing,
            bVsites,bVsiteAromatics,forcefield,ffdir,
            mHmult,nssbonds,ssbonds,
            long_bond_dist,short_bond_dist,bDeuterate,bChargeGroups,bCmap,
            bRenumRes,bRTPresname);
    
    if (!cc->bAllWat)
      write_posres(posre_fn,pdba,posre_fc);

    if (bITP)
      gmx_fio_fclose(itp_file);

    /* pdba and natom have been reassigned somewhere so: */
    cc->pdba = pdba;
    cc->x = x;
    
    if (debug) {
      if (cc->chainid == ' ')
        sprintf(fn,"chain.pdb");
      else
        sprintf(fn,"chain_%c.pdb",cc->chainid);
      cool_quote(quote,255,NULL);
      write_sto_conf(fn,quote,pdba,x,NULL,ePBC,box);
    }
  }

  if (watermodel == NULL) {
    for(chain=0; chain<nch; chain++) {
      if (chains[chain].bAllWat) {
        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.");
      }
    }
  } else {
    sprintf(buf_fn,"%s%c%s.itp",ffdir,DIR_SEPARATOR,watermodel);
    if (!fflib_fexist(buf_fn)) {
      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.",
                buf_fn,watermodel);
    }
  }

  print_top_mols(top_file,title,ffdir,watermodel,nincl,incls,nmol,mols);
  gmx_fio_fclose(top_file);

  gmx_residuetype_destroy(rt);
    
  /* now merge all chains back together */
  natom=0;
  nres=0;
  for (i=0; (i<nch); i++) {
    natom+=chains[i].pdba->nr;
    nres+=chains[i].pdba->nres;
  }
  snew(atoms,1);
  init_t_atoms(atoms,natom,FALSE);
  for(i=0; i < atoms->nres; i++)
    sfree(atoms->resinfo[i].name);
  sfree(atoms->resinfo);
  atoms->nres=nres;
  snew(atoms->resinfo,nres);
  snew(x,natom);
  k=0;
  l=0;
  for (i=0; (i<nch); i++) {
    if (nch>1)
      printf("Including chain %d in system: %d atoms %d residues\n",
             i+1,chains[i].pdba->nr,chains[i].pdba->nres);
    for (j=0; (j<chains[i].pdba->nr); j++) {
      atoms->atom[k]=chains[i].pdba->atom[j];
      atoms->atom[k].resind += l; /* l is processed nr of residues */
      atoms->atomname[k]=chains[i].pdba->atomname[j];
      atoms->resinfo[atoms->atom[k].resind].chainid = chains[i].chainid;
      copy_rvec(chains[i].x[j],x[k]);
      k++;
    }
    for (j=0; (j<chains[i].pdba->nres); j++) {
      atoms->resinfo[l] = chains[i].pdba->resinfo[j];
      if (bRTPresname) {
        atoms->resinfo[l].name = atoms->resinfo[l].rtp;
      }
      l++;
    }
  }
  
  if (nch>1) {
    fprintf(stderr,"Now there are %d atoms and %d residues\n",k,l);
    print_sums(atoms, TRUE);
  }
  
  fprintf(stderr,"\nWriting coordinate file...\n");
  clear_rvec(box_space);
  if (box[0][0] == 0) 
    gen_box(0,atoms->nr,x,box,box_space,FALSE);
  write_sto_conf(ftp2fn(efSTO,NFILE,fnm),title,atoms,x,NULL,ePBC,box);

  printf("\t\t--------- PLEASE NOTE ------------\n");
  printf("You have successfully generated a topology from: %s.\n",
         opt2fn("-f",NFILE,fnm));
  if (watermodel != NULL) {
    printf("The %s force field and the %s water model are used.\n",
           ffname,watermodel);
  } else {
    printf("The %s force field is used.\n",
           ffname);
  }
  printf("\t\t--------- ETON ESAELP ------------\n");
  

  thanx(stdout);
  
  return 0;
}