# include "../src-lib/dist.h"
# include "traj.h"

# define MAX_MOLS 100


/*#####################################################################################*/
/*                          Stuff for get coord mdb                                    */
/*#####################################################################################*/
# define BiG_NuM_MaX 30000
# define INTSIZE 4
# define EOR -2                    /* defines End Of Record */

static int Recsize, Recpos;
static short Ixk[BiG_NuM_MaX];
static float Coords[BiG_NuM_MaX], Sclp, Sclw, Sclp1;
/*#####################################################################################*/



/*=======================================================================================
                Cube Subroutines for use with Encad Trajectories
=======================================================================================*/


/* got this from britt */
AtomList **** alloc3dAtomListptr (int x, int y, int z)

{
  int i,j;
  AtomList ****ret;
  AtomList ***array;
  char *buff;
  char *data;
  int size,memneed;
  memneed=y*sizeof(AtomList **)+z*y*sizeof(AtomList *);
  size=memneed*x+x*sizeof(AtomList ***);
  data=(char *)calloc(size,1);
  ret=(AtomList ****)data;
  for (j=0; j<x; j++)
  {
    buff=&data[x*sizeof(AtomList ***)+memneed*j];
    array=(AtomList ***)buff;
    buff+=y*sizeof(AtomList **);
    for (i=0; i<y; i++)
    {
      array[i]=(AtomList **)buff;
      buff+=z*sizeof(AtomList *);
    }
    ret[j]=array;
  }
  return ret;
}


/* returns pointers to Unit_Cube in three dimensions x, y, & z */
/*        ============                                                      */
theCube * SetuptheCube (file_records * f, float * bxsidel, double nb_cutoff)
/*        ============                                                      */
{
  int ii;

  theCube * C = (theCube *)calloc(1,sizeof(theCube));

  /* non-bonded cutoff */
  C->nb_cut = nb_cutoff;

  /* defining the number of bins */
  C->nbinx = floor(bxsidel[0]/nb_cutoff);
  C->nbiny = floor(bxsidel[1]/nb_cutoff);
  C->nbinz = floor(bxsidel[2]/nb_cutoff);

  /* defining x, y, z size of each unit box           */
  /* using definition defined for minimum coordinates */
  C->sbinx = bxsidel[0]/C->nbinx;
  C->sbiny = bxsidel[1]/C->nbiny;
  C->sbinz = bxsidel[2]/C->nbinz;

  /* print out cube divisions and # of divisions */
  fprintf(stderr, "# non bonded cutoff for unit boxes is %.2lf\n", nb_cutoff);
  fprintf(stderr, "# for x dimension: %10.5lf gives %d bins of %10.5lf A.\n",
	 bxsidel[0], C->nbinx, C->sbinx);
  fprintf(stderr, "# for y dimension: %10.5lf gives %d bins of %10.5lf A.\n",
	 bxsidel[1], C->nbiny, C->sbiny);
  fprintf(stderr, "# for z dimension: %10.5lf gives %d bins of %10.5lf A.\n", 
	 bxsidel[2], C->nbinz, C->sbinz);
  
  /* using number of bins, make 3d array of AtomList pointers */
  C->BinnedAtoms = alloc3dAtomListptr(C->nbinx, C->nbiny, C->nbinz);

  /* for each atom in pdb, make AtomList */
  C->theatoms = (AtomList *)calloc(f->atomnum, sizeof(AtomList));
  for(ii=0;ii<f->atomnum;ii++)  {
#ifdef OFF
    atom_record * a = f->atoms[ii];
    if(!strcmp("HOH", a->resname))
      fprintf(stderr, "%s %s %d %f %f %f\n",  
	      a->resname, a->name, a->number,a->x, a->y, a->z);
#endif
 
    C->theatoms[ii].atom = f->atoms[ii];    
  }

  C->atomnum = f->atomnum;

  return C;
}


/* these two should be called for each trajectory analyzed */
void BinMoleculesIntheCube (file_records * f, theCube * C, float * box_min) 
{
  int 
    ii, x_reg, y_reg, z_reg;

  AtomList
    * findend, * new,
    **** A = C->BinnedAtoms;

  for(ii=0; ii<f->atomnum; ii++)  {
    atom_record
      * b,
      * a = f->atoms[ii];

    x_reg = floor((a->x - box_min[0])/C->sbinx);
    y_reg = floor((a->y - box_min[1])/C->sbiny);
    z_reg = floor((a->z - box_min[2])/C->sbinz);

#ifdef OFF
    printf("-- %d %d %d %d\n", x_reg, y_reg, z_reg, f->atomnum);
    printf("-- %lf %lf %lf\n", a->x - box_min[0], a->y - box_min[1], a->z - box_min[2]);
#endif

    /* must find a better way to do */
    if(A[x_reg][y_reg][z_reg] == NULL) {
      A[x_reg][y_reg][z_reg] = &(C->theatoms[ii]);

#ifdef OFF
      b = A[x_reg][y_reg][z_reg]->atom;
      printf("#1 %s %s %d %d\n", b->name, b->resname, b->number, b->resnum);
#endif
    }
    else {
      findend = A[x_reg][y_reg][z_reg];
      while( findend->NextAtom != NULL)  {

#ifdef OFF
	b = findend->atom;
	printf("%s %s %d %d\n", b->name, b->resname, b->number, b->resnum);
#endif
	findend = findend->NextAtom;

      }
      findend->NextAtom = &(C->theatoms[ii]);

#ifdef OFF
      b = findend->NextAtom->atom;
      printf("%s %s %d %d\n", b->name, b->resname, b->number, b->resnum);
#endif
    }    
  }
}


/* this routine resets all the AtomList pointers in theCube to NULL */
void ResettheCube (theCube * C)
{
  int
    xx, yy, zz;

  for(xx=0;xx<C->nbinx;xx++)  {
    for(yy=0;yy<C->nbiny;yy++)  {
      for(zz=0;zz<C->nbinz;zz++)  {
	C->BinnedAtoms[xx][yy][zz] = NULL;
      }
    }
  }

  for(xx=0;xx<C->atomnum;xx++)  {
    C->theatoms[xx].NextAtom = NULL;
  }
}



/* something is fishy about this */
void FindNeighbors (file_records *f, theCube * C, float * box_min)
{
  int 
    ii, 
    xreg, yreg, zreg,
    xx, yy, zz,
    startx, starty, startz,
    stopx, stopy, stopz,
    nx, ny, nz;

  AtomList * thit;

  startx = C->nbinx-1;
  starty = C->nbiny-1;
  startz = C->nbinz-1;

  stopx = C->nbinx+2;
  stopy = C->nbiny+2;
  stopz = C->nbinz+2;
 
  for(ii=0;ii<f->atomnum && ii<1 ;ii++)  {
    atom_record * a = f->atoms[ii];

    printf("for %s %s %d %d\n", a->name, a->resname, a->number, a->resnum);


    xreg = (int)floor((a->x - box_min[0])/C->sbinx);
    yreg = (int)floor((a->y - box_min[1])/C->sbiny);
    zreg = (int)floor((a->z - box_min[2])/C->sbinz);

    printf("%d %d %d %d *\n",  ii, xreg, yreg, zreg);

    for(xx=startx; xx<stopx; xx++)  {
      for(yy=starty; yy<stopy; yy++)  {
	for(zz=startz; zz<stopz; zz++)  {
	  nx = (xreg + xx)%C->nbinx;
	  ny = (yreg + yy)%C->nbiny;
	  nz = (zreg + zz)%C->nbinz;

	  printf("%d %d %d %d\n",  ii, nx, ny, nz);

	  thit = C->BinnedAtoms[nx][ny][nz];
	  while(thit!= NULL)  {
	    atom_record * b = thit->atom;
	    printf("%s %s %d %d\n", b->name, b->resname, b->number, b->resnum);
	    thit = thit->NextAtom;
	  }

	  if(nx<0 || nx >= C->nbinx ||
	     ny<0 || ny >= C->nbiny ||
	     nz<0 || nz >= C->nbinz)
	    puts("Got one!!");
	}
      }
    }
  }
}

/* ====================== End of Cube Routines =========================== */




/*     A General Histogram Program       */
/*   ==============                      */
void BasicHistogram 
/*   ==============                      */
(double *input, int numinput, double min, double max, double binwidth, int numofbins)
{
  int ii, jj;

  double * histogram;

  if (min==max)  {
    min = max = input[0];

    for(ii=1; ii<numinput; ii++)  {

      if (input[ii] < min) {
	min = input[ii];
	continue;
      }

      if (input[ii] > max) {
	max = input[ii];
	continue;
      }
    }
    printf("# found min %5.3le and max %5.3le\n", min, max);
  }
  else
    printf("# min %6.3lf ; max %6.3lf\n", min, max);


  if(numofbins == 0 && binwidth != 0)  {
    numofbins = ceil((max-min)/binwidth);
    printf("# binwidth %.2le defines %d bins\n", binwidth, numofbins);
  }
  else if (numofbins != 0 && binwidth == 0)  {
    binwidth = (1.0001*max-min)/numofbins;
    printf("# %d bins yields %6.3lf binwidth\n", numofbins, binwidth);
  }
  else if (numofbins == 0 && binwidth == 0)  {
    numofbins = ceil(sqrt(numinput));
    binwidth = (1.0001*max-min)/numofbins;
    printf("# used sqrt of number of input values to find %d bins and %6.3le binwidth\n", 
	   numofbins, binwidth);
  }
  else
    printf("# using values %d bins and %6.3lf binwidth\n", numofbins, binwidth);


  histogram = (double *)calloc(numofbins, sizeof(double));

  for (ii=0; ii<numinput; ii++)  {
    jj = floor((input[ii] - min)/binwidth);
    histogram[jj]++;
  }

  printf("# bins            values\n");
  printf("# -----------     ---------\n");
  for (ii=0; ii<numofbins; ii++)
    printf("  %.3le           %9.2lf\n", ii+1, ii*binwidth, histogram[ii]);

}




/*   ================                                                 */
void VaryBySideLength (double num, double * oto, float bxsidel)
/*   ================                                                 */
{
  int ii, jj;

  for(ii=-1, jj=0; jj<3; ii++, jj++)  {
    oto[jj] = num + (ii*bxsidel);
  }
}


/*---------------------------------------------------------------------------*/
/*         calculates 27 possible distances between two atoms                */
/*       one centered in the box; the other reflected 27 times               */
/*---------------------------------------------------------------------------*/
/*       ============                                                        */
double * theDistances
/*       ============                                                        */
(atom_record * center, atom_record * radial, float * bxsidel)
{
  int
    ii,
    jj,
    kk,
    count = 0;

  double
    xx[3], yy[3], zz[3],
    distance[27],
    dx = radial->x - center->x,
    dy = radial->y - center->y,
    dz = radial->z - center->z;

  dx = dx - (bxsidel[0] * CROUND(dx/bxsidel[0]));
  dy = dy - (bxsidel[1] * CROUND(dy/bxsidel[1]));
  dz = dz - (bxsidel[2] * CROUND(dz/bxsidel[2]));
  
  VaryBySideLength (dx, xx, bxsidel[0]);
  VaryBySideLength (dy, yy, bxsidel[1]);
  VaryBySideLength (dz, zz, bxsidel[2]);

  for(ii=0; ii<3; ii++)  {
    for(jj=0; jj<3; jj++)  {
      for(kk=0; kk<3; kk++)  {
	distance[count] = sqrt(xx[ii]*xx[ii] + yy[jj]*yy[jj] + zz[kk]*zz[kk]);
	count++;
      }
    }
  }
  return distance;
}


/* this keeps the integrity of entire molecule for boxing  */
/*   ====================                                  */
void BoxItWholeMolecules (file_records * f, float *bxsidel)
/*   ====================                                  */
{
  int ii, resnum=-1;
  float dx, dy, dz;
  char *resname = "new";

  for(ii=0; ii<f->atomnum; ii++)  {
    atom_record * a = f->atoms[ii];

    if(resnum != a->resnum || strcmp(resname, a->resname))  {
      dx = (bxsidel[0] * CROUND(a->x/bxsidel[0]));
      dy = (bxsidel[1] * CROUND(a->y/bxsidel[1]));
      dz = (bxsidel[2] * CROUND(a->z/bxsidel[2]));
      resnum = a->resnum; strcpy(resname, a->resname);
    }

    a->x -= dx; a->y -= dy; a->z -= dz;
  }
}



/* try to build in asymmetric min and max for sides */
/*   =====                                                 */
void BoxIt (file_records * f, float *bxsidel)
/*   =====                                                 */
{
  int ii;

  for(ii=0; ii<f->atomnum; ii++)  {
    atom_record * a = f->atoms[ii];

    /* ie. boxside = 18, -10 -->  -10 - (18*CROUND(-10/18)) */
    /*        new coordinate is -10 - (18*-1) = 8           */
    a->x = a->x - (bxsidel[0] * CROUND(a->x/bxsidel[0]));
    a->y = a->y - (bxsidel[1] * CROUND(a->y/bxsidel[1]));
    a->z = a->z - (bxsidel[2] * CROUND(a->z/bxsidel[2]));


    /*----------------------------------------------------
        A macro from mark, look for it in a header file
      
               CROUND(x) = ANINT(x) for all x
             2.4 becomes  2      2.5 becomes  3
            -2.4 becomes -2     -2.5 becomes -3
      ----------------------------------------------------*/
  }
}


/*   ===============                                                */
void FindSideLengths (float * min, float * max, float * bxsidel)
/*   ===============                                                */
{

  bxsidel[0] = max[0] - min[0];
  bxsidel[1] = max[1] - min[1];
  bxsidel[2] = max[2] - min[2];

  fprintf(stderr, "# used following dimensions for the box.\n");
  fprintf(stderr, "#\t\t  x\t  y\t  z\n");
  fprintf(stderr, "# min\t\t%6.2f\t%6.2f\t%6.2f\n", min[0], min[1], min[2]);
  fprintf(stderr, "# max\t\t%6.2f\t%6.2f\t%6.2f\n", max[0], max[1], max[2]);
  fprintf(stderr, "# bxsidel\t%6.2f\t%6.2f\t%6.2f\n#\n", bxsidel[0], bxsidel[1], bxsidel[2]);

}


/*     ===============                                                      */
double FindMinDistance 
/*     ===============                                                      */
(atom_record * a, atom_record * b, float * bxsidel)
{
  double 
    r2dist, 
    dx = a->x - b->x, 
    dy = a->y - b->y,
    dz = a->z - b->z;

  /*       for minimum image convention       */
  dx = dx - (bxsidel[0]*CROUND(dx/bxsidel[0]));
  dy = dy - (bxsidel[1]*CROUND(dy/bxsidel[1]));
  dz = dz - (bxsidel[2]*CROUND(dz/bxsidel[2]));

  r2dist = dx*dx + dy*dy + dz*dz; 

  return sqrt(r2dist);
}



/*###################################################################################
  This is the new trajectory reading program,  for a test see read_mdb.sgi
  in ~/code/read_mdb/ .
  incorporated britt's code
###################################################################################*/


/* Gets the int size of a byte record from record header */
/* returns EOF if at End Of File                         */
static int rget_int(FILE *mdb, int *val)
{
  int i, t;
  char *c;
  c=(char *)val;
  for(i=0;i<INTSIZE;i++)  {
    t=getc(mdb);
    if(t==EOF) 
      return EOF;
    c[i]=t;
  }
  return 0;
}


/* stores the size of the record in Recsize and sets Recpos to zero */
/* returns EOF if at End Of File                                    */
int newrecord(FILE *mdb)
{
  int stat;
  stat=rget_int(mdb, &Recsize);           /* read the record header */
  Recpos=0;
  return stat;
}


/* goes to the end of a record */
void endrecord(FILE *mdb)
{
  int t;
  while (Recpos<Recsize)  {
    t=getc(mdb);
    Recpos++;
  }
  rget_int(mdb, &t);                       /* what is this for??? */
}


/* reads a byte */
int ftget_byte(FILE *mdb, char *val)
{
  int t;
  t=getc(mdb);
  if (t==EOF)
    return EOF;
  Recpos++;
  val[0]=t;
  return 0;
}


/* reads a binary integer */
int ftget_int(FILE *mdb, int *val)
{
  int i, t;
  char *c;
  c=(char *)val;
  for(i=0;i<INTSIZE;i++)  {
    if(Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if(t==EOF)
      return EOF;
    c[i]=t;
    Recpos++;
  }
  return 0;
}


/* read in a short integer */
int ftget_short(FILE *mdb, short *val)
{
  int i,t;
  char *c;
  c=(char *)val;
  for (i=0; i<2; i++)  {
    if (Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if (t==EOF)
      return EOF;
    c[i]=t;
    Recpos++;
  }
  return 0;
}


/* read in a float */
int ftget_float(FILE *mdb, float *val)
{
  int i,t;
  char *c;
  c=(char *)val;
  for (i=0; i<4; i++)  {
    if (Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if (t==EOF)
      return EOF;
    c[i]=t;
    Recpos++;
  }
  return 0;
}


/* do i need this or not? */
/* read a string */
int ftget_str(FILE *mdb, int n, char *p)
{
  int i,t;
  for (i=0; i<n; i++)  {
    if (Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if (t==EOF)
      return EOF;
    p[i]=t;
    Recpos++;
  }
  p[i]==0;       /* Cap it off */
  return 0;
}



/* routine reads in a mdb file     */
/* returns EOF at the End Of File  */
/* otherwise returns 0             */
/*  ==============                                                          */
int get_mdb_coords (file_records *f, int * istep, float * min, float * max)
/*  ==============                                                          */

{
  int icode, nat3p, n, i, j;
  static int after_first_record = 0;
  signed char ixb;
  static FILE * mdb;


  if(after_first_record)  {

    if (newrecord(mdb)==EOF)
      return EOF;
    ftget_int(mdb,&icode);      /* icode specifies stored coords' format */
    endrecord(mdb);
    
    newrecord(mdb);
    ftget_int(mdb,istep);
    ftget_int(mdb,&n);
    if (icode==0) { /* icode = 0 indicates trajectory stored as a byte shift */
      for (i=0; i<n; i++)  {
	ftget_byte(mdb,(char *)&ixb);
	Ixk[i]+=ixb;
      }
    }
    else  {        /* otherwise coordinates written out completely in a short */
      for (i=0; i<n; i++)
	ftget_short(mdb,&(Ixk[i]));
    }
    endrecord(mdb);
    
    for (i=0; i<n; i+=3)  {       /* expands coordinates into floats */
      atom_record * a = f->atoms[i/3];

      a->x = (double)(Ixk[i])*Sclp1*0.01;
      a->y = (double)(Ixk[i+1])*Sclp1*0.01;
      a->z = (double)(Ixk[i+2])*Sclp1*0.01;
    }
  }
  /* this reads in first record in mdb */
  /* sets Sclp and Sclp1               */
  /* gets min and max arrays           */
  else {  

    if(!(mdb=fopen("fort.1", "r"))) {
      STDERR("Couldn't open fort.1");
      exit(1);
    }

    newrecord(mdb);

    ftget_int(mdb,istep);
    ftget_int(mdb,&n);

    for (i=0; i<n; i++)       /* getting raw coordinates */
      ftget_short(mdb,&(Ixk[i]));

    ftget_int(mdb,&nat3p);     /* old vs new style coords (see getix2) */
    ftget_float(mdb,&Sclp);    /* scale to multiply coordinates        */
    Sclp1=1.0/Sclp;           
    ftget_float(mdb,&Sclw);    /* something not used                   */

    /* getting min and max for x, y, & z */
    for (i=0; i<3; i++)
      ftget_float(mdb, &min[i]);
    for (i=0; i<3; i++) 
      ftget_float(mdb, &max[i]);

    /* expands coordinates into floats and writes into file record */
     for (i=0; i<n; i+=3)  {       /* expands coordinates into floats */
      atom_record * a = f->atoms[i/3];

      a->x = (double)(Ixk[i])*Sclp1*0.01;
      a->y = (double)(Ixk[i+1])*Sclp1*0.01;
      a->z = (double)(Ixk[i+2])*Sclp1*0.01;
    }
    endrecord(mdb);   

    after_first_record++;
  }  

  return 0;
}


#ifdef OFF
/*********************************************************************************/
/* use the above instead , get_mdb_coords */
/*********************************************************************************/


/*-----------------------------------------------------------------------------
 *  uses Mike's program getix2
 * 
 * This function returns:
 *
 *         a. template file_record with new xyz coordinates at a particular
 *            istep
 *         b. integer which indicates if reached end of trajectory file:
 *            if integer <= 0. in getix2.lc integer is -1 to indicate end
 *            of trajectory file.
 *         c. time in steps, where 1 step is 2 femtoseconds.  
 *         d. length of half the side of the box (should check if positive)
--------------------------------------------------------------------------------*/

/*   ====================                                    */
long read_coord_from_traj
/*   ====================                                    */
(file_records * f, long * time, float * min, float * max)
{
  
  short ix[30000];
  char  ixb[30000];

  float
    xmax[3],
    xmin[3];
  
  long
    iunit = 1,
    istep,
    n,
    i,
    j;
  
  
  /* Mike's prgram that reads in trjectory files
   *    - requires that input file be named fort.1
   *    - has internal marker so that each time it is called it outputs next
   *      istep
   */
  
  getix2_ ( &iunit, &istep, &n, ix, ixb, xmin, xmax );

  *time = istep;

  min[0] = xmin[0];
  min[1] = xmin[1];
  min[2] = xmin[2];

  max[0] = xmax[0];
  max[1] = xmax[1];
  max[2] = xmax[2];
  
  if ( n <= 0 )  {
    printf("# ");
    STDOUT("The trajectory is exhausted.");
    return n ;
  }

  else {
    for ( i=0; i < n; i+=3 ) {
      atom_record * a = f->atoms[i/3];

      a->x = 0.01*ix[i];
      a->y = 0.01*ix[i+1];
      a->z = 0.01*ix[i+2];
    }

    return n;
  }
}

#endif



/*   =========                                                           */
void PrintPtrs (atom_record ** list, int num)
/*   =========                                                           */
{
  file_records 
    hoh, 
    *oho = &hoh;

  oho->atomnum = num;
  oho->atoms = list;

  fprintf(stderr, "# for a set of pointers: \n");
  write_pdb_file (stderr, oho, 0);
}


/*   ==========                                                           */
void AssignPtrs (file_records * f, atom_record ** list, Ptr_Info * p)
/*   ==========                                                           */
{
  int
    count,
    num = 0;

  for (count=0; count<f->atomnum; count++)  {
    atom_record * a = f->atoms[count];
    
    if( !strcmp(a->resname, p->pdb_nm)  && 
       (!strcmp(a->name, p->atm1) || !strcmp(a->name, p->atm2)) )  {
      list[num] = a;
      num++;
    }
  }
}


/*  ==========                                                    */
int count_atom (file_records * f, char * resname, char * name)
/*  ==========                                                    */
{
  int
    num = 0,
    count;
  
  for (count=0; count<f->atomnum; count++)  {
    atom_record * a = f->atoms[count];
    
    if ( !strcmp(a->resname, resname)  && !strcmp(a->name, name) )
      num++;
  }

  if (num == 0)
    fprintf(stderr, "\n\nHey Buckaroo! don't have any %s in the file!\n", resname);

  return num;
}


int FindMoleculesInPDB (file_records * f, char ** mols)
{
  int ii, jj=0, kk, new_molecule;

  mols[jj] = f->atoms[jj]->resname;

  for(ii=1; ii<f->atomnum; ii++)  {
    char * a = f->atoms[ii]->resname;

    for (kk=0, new_molecule=1; kk<jj+1; kk++)  {
      if ( !strcmp(a, mols[kk]) ) {
	new_molecule=0;
	break;
      }
    }

    if ( new_molecule )  {
      jj++;
      mols[jj]=a;
    }
    
    if (jj>MAX_MOLS)  {
      fprintf(stderr, "Too Many Different Molecules!!\n");
      fprintf(stderr, "Will print out what i have found\n");
      break;
    }
  }
  
  return ++jj;
}



int CheckMolecule (file_records * f, char * resname) 
{
  int count;

  for (count=0; count<f->atomnum; count++)  {
    atom_record * a = f->atoms[count];
    
    if ( !strcmp(a->resname, resname) )
      return 1;
  }
  
  fprintf(stderr, "[Hey Buckaroo! don't have any %s in the file!]\n", resname);
  printf("# Hey Buckaroo! didn't find any %s in the file!\n", resname);
  return 0;
}


/*         ==================                                                      */
Ptr_Info **CreatePtrInfoArray (int argc, char **argv, file_records * f, int *num_A)
/*         ==================                                                      */
{
  int 
    jj,
    ii,
    count,
    num,
    num_molecules;

  char 
    * hoh = "HOH",
    * ure = "URE",
    * act = "ACT",
    * prp = "PRP",
    * poh = "POH",
    * atm = "ATM",
    * ibe = "IBE",
    * bnz = "BNZ",
    * but = "BUT",
    * guh = "GUH",
    * cl  = "CL ",
    * mth = "MTH",
    * molecules[MAX_MOLS];

  Ptr_Info ** Array;

  num = 0;
  if(ArgvMatch(argc, argv, "hoh") && CheckMolecule(f, "HOH"))  {
    molecules[num] = hoh;
    num++;
  }
  if(ArgvMatch(argc, argv, "ure") && CheckMolecule(f, "URE"))  {
    molecules[num] = ure;
    num++;
  }
  if(ArgvMatch(argc, argv, "act") && CheckMolecule(f, "ACT"))  {
    molecules[num] = act;
    num++;
  }
  if(ArgvMatch(argc, argv, "prp") && CheckMolecule(f, "PRP"))  {
    molecules[num] = prp;
    num++;
  }
  if(ArgvMatch(argc, argv, "poh") && CheckMolecule(f, "POH"))  {
    molecules[num] = poh;
    num++;
  }
  if(ArgvMatch(argc, argv, "atm") && CheckMolecule(f, "ATM"))  {
    molecules[num] = atm;
    num++;
  }
  if(ArgvMatch(argc, argv, "ibe") && CheckMolecule(f, "IBE"))  {
    molecules[num] = ibe;
    num++;
  }
  if(ArgvMatch(argc, argv, "bnz") && CheckMolecule(f, "BNZ"))  {
    molecules[num] = bnz;
    num++;
  }
  if(ArgvMatch(argc, argv, "but") && CheckMolecule(f, "BUT"))  {
    molecules[num] = but;
    num++;
  }
  if(ArgvMatch(argc, argv, "guh") && CheckMolecule(f, "GUH"))  {
    molecules[num] = guh;
    num++;
  }
  if(ArgvMatch(argc, argv, "cl ") && CheckMolecule(f, "CL "))  {
    molecules[num] = cl;
    num++;
  }
  if(ArgvMatch(argc, argv, "mth") && CheckMolecule(f, "MTH"))  {
    molecules[num] = mth;
    num++;
  }
  
  if (num == 0)  {

    num_molecules = FindMoleculesInPDB(f, molecules);

    if (num_molecules > 0)
      fprintf(stderr, "# Have found %d types of molecules in the input pdb:\n",
	      num_molecules);
    else
      fprintf(stderr, "# Have found only one type of molecule in the input pdb:\n");

    for(ii=0; ii<num_molecules; ii++)
      fprintf(stderr, "# \t%s\n", molecules[ii]);

    num = num_molecules;
  }

  if(num==0)  {
    fprintf(stderr, "didn't find any molecules to work on!!\n");
    STDERR("\n\nPlease specify molecules over which to calculate RDF.");
    fprintf(stderr, "Molecules so far in database:\n\n");
    fprintf(stderr, "\t-hoh for water\n");
    fprintf(stderr, "\t-ure for urea\n");
    fprintf(stderr, "\t-act for acetone\n");
    fprintf(stderr, "\t-poh for isopropanol\n");
    fprintf(stderr, "\t-prp for propane\n");
    fprintf(stderr, "\t-atm for acetamide\n");
    fprintf(stderr, "\t-ibe for isobutylene\n");
    fprintf(stderr, "\t-bnz for benzene\n");
    fprintf(stderr, "\t-but for butane\n");
    fprintf(stderr, "\t-guh for guanidine\n");
    fprintf(stderr, "\t-cl for chlorine\n");
    fprintf(stderr, "\t-mth for methane\n");
    fprintf(stderr, "\n");
    exit(-1);
  }

/*  fprintf(stderr, "# num = %d\n", num); */
  
  Array = (Ptr_Info **)calloc(num, sizeof(Ptr_Info *)); 
  for(count=0; count<num; count++) 
    Array[count] = (Ptr_Info *)malloc(sizeof(Ptr_Info));


  for(count=0, ii=0; count<num; count++)  {
    if( !strcmp(molecules[count], hoh) ) {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "water");
      strcpy(p->pdb_nm, "HOH");
      strcpy(p->atm1, "OH ");
      strcpy(p->atm2, "zzz");
      strcpy(p->ctr_atm, "OH ");
      p->mod = 1;
      p->atm_nm = 3;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], ure) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "urea");
      strcpy(p->pdb_nm, "URE");
      strcpy(p->atm1, "N1 ");
      strcpy(p->atm2, "N6 ");
      strcpy(p->ctr_atm, "C4 ");
      p->mod = 2;
      p->atm_nm = 8;
      ii++;
      continue;
    }	
    if( !strcmp(molecules[count], act) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "acetone");
      strcpy(p->pdb_nm, "ACT");
      strcpy(p->atm1, "C1 ");
      strcpy(p->atm2, "C7 ");
      strcpy(p->ctr_atm, "C5 ");
      p->mod = 2;
      p->atm_nm = 10;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], prp) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "propane");
      strcpy(p->pdb_nm, "PRP");
      strcpy(p->atm1, "CA ");
      strcpy(p->atm2, "CC ");
      strcpy(p->ctr_atm, "CB ");
      p->mod = 2;
      p->atm_nm = 8;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], poh) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "isopropanol");
      strcpy(p->pdb_nm, "POH");
      strcpy(p->atm1, "C1 ");
      strcpy(p->atm2, "C9 ");
      strcpy(p->ctr_atm, "C5 ");
      p->mod = 2;
      p->atm_nm = 12;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], atm) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "acetamide");
      strcpy(p->pdb_nm, "ATM");
      strcpy(p->atm1, "C1 ");
      strcpy(p->atm2, "N7 ");
      strcpy(p->ctr_atm, "C5 ");
      p->mod = 2;
      p->atm_nm = 9;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], ibe) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "isobutylene");
      strcpy(p->pdb_nm, "IBE");
      strcpy(p->atm1, "C1 ");
      strcpy(p->atm2, "C9 ");
      strcpy(p->ctr_atm, "C5 ");
      p->mod = 2;
      p->atm_nm = 12;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], guh) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "guanidine");
      strcpy(p->pdb_nm, "GUH");
      strcpy(p->atm1, "N1 ");
      strcpy(p->atm2, "N8 ");
      strcpy(p->ctr_atm, "C4 ");
      p->mod = 2;
      p->atm_nm = 10;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], cl) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "chlorine");
      strcpy(p->pdb_nm, "CL ");
      strcpy(p->atm1, "CL ");
      strcpy(p->atm2, "zzz");
      strcpy(p->ctr_atm, "CL ");
      p->mod = 1;
      p->atm_nm = 1;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], mth) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "methane");
      strcpy(p->pdb_nm, "MTH");
      strcpy(p->atm1, "CA ");
      strcpy(p->atm2, "CA ");
      strcpy(p->ctr_atm, "CA ");
      p->mod = 1;
      p->atm_nm = 5;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], bnz) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "benzene");
      strcpy(p->pdb_nm, "BNZ");
      strcpy(p->atm1, "C1 ");
      strcpy(p->atm2, "C2 ");
      strcpy(p->ctr_atm, "C3 ");
      p->mod = 2;
      p->atm_nm = 12;
      ii++;
      continue;
    }
    if( !strcmp(molecules[count], but) )  {
      Ptr_Info * p = Array[ii];
      strcpy(p->name, "butane");
      strcpy(p->pdb_nm, "BUT");
      strcpy(p->atm1, "CA ");
      strcpy(p->atm2, "CB ");
      strcpy(p->ctr_atm, "CC ");
      p->mod = 2;
      p->atm_nm = 14;
      ii++;
      continue;
    }

    printf("# did not have record of %s in the database\n", molecules[count]);
  } 

  if( ii==0 )  {
    fprintf(stderr, "did not recognize any molecules in database!!");
    exit(-2);
  }

/*  fprintf(stderr, "# ii = %d\n", ii); */
  *num_A = ii;
  return Array;
}


/*-----------------------------------------------------------------

  function formats mike's pdb file into more regular pdb format
  what it does:
     1. for waters:
          - residue name changes from W10 to HOH
	  - residue numbers increase every water molecule instead 
	    of every 10 molecules
	  - formats all oxygens to OW
	  - formats all hydrogens to H1 or H2 depending on position


     2. for ureas:
          - formats residue numbers so that they are sequential

     3. the END of the protein gets the same resname as previous atom.
        also, for it is the c termini's oxygen, OXT instead of OT.
	and same residue number as the previous atom.

     4. flags all the hydrogens as deselected

     5. only dealt with cl ion and na ion

     6. added same ability to handle acetones as for ureas

------------------------------------------------------------------*/
/*   ==========                                                   */
void pdb_format (file_records * f)
/*   ==========                                                   */
{
  int
    atmcount, 
    protrescnt        = 0, pres = -1,
    water_count       = 1, water_num         = 1,
    urea_count        = 1, urea_num          = 1,
    acetone_count     = 1, acetone_num       = 1,
    propane_count     = 1, propane_num       = 1,
    isopropanol_num   = 1, isopropanol_count = 1,
    acetamide_num     = 1, acetamide_count   = 1,
    isobutylene_num   = 1, isobutylene_count = 1,
    guh_num           = 1, guh_count         = 1,
    mth_num           = 1, mth_count         = 1,
    bnz_num           = 1, bnz_count         = 1,
    but_num           = 1, but_count         = 1,
    cl_count_num      = 1, na_count_num      = 1;
  

  for (atmcount=0; atmcount<f->atomnum; atmcount++)  {
    atom_record * a = f->atoms[atmcount];

    /* insuring protein residue number is correct */
    if ( PROTEIN_ATOM_P(a) && strcmp(a->resname, "END") ) {
      if(pres != a->resnum)  {protrescnt++; pres = a->resnum;}
      a->resnum = protrescnt;
      continue;
    }

    else if(!strcmp(a->resname, "END")) {
      strcpy(a->resname, f->atoms[atmcount-1]->resname);
      a->resnum = f->atoms[atmcount-1]->resnum;
    }

    /* addressing the waters */
    else if (!strcmp(a->resname,"WAT") || !strcmp(a->resname, "HOH") 
	     || !strcmp(a->resname,"DOD") || !strcmp(a->resname,"WA1")
	     || !strcmp(a->resname,"W10")) {
      strcpy(a->resname, "HOH");
      a->resnum = water_num;

      if (a->name[0] == 'O')   
	strncpy(a->name, "OH ", 4);
      
      
      if (a->name[0] == 'H')  {
	if (f->atoms[atmcount-1]->name[0] == 'O')
	  strncpy(a->name, "H1 ", 4);
	else if (f->atoms[atmcount-1]->name[0] == 'H')
	  strncpy(a->name, "H2 ", 4);
	
	a->extraflag = ' ';
      }
      
      if ( (water_count%3) == 0 )
	water_num++;

      water_count++;
    }
      
    /* making sure resnum are sequential for urea */
    else if ( !strcmp(a->resname, "URE") )  {
      a->resnum = urea_num;
      
      if ( (urea_count%8) == 0 )
	urea_num++;
      
      urea_count++;
    }
      
    /* making sure resnum are sequential for acetones */
    else if ( !strcmp(a->resname, "ACT") )  {
      a->resnum = acetone_num;
      
      if ( (acetone_count%10) == 0 )
	acetone_num++;
      
      acetone_count++;
    }
      
    /* making sure resnum are sequential for propanes */
    else if ( !strcmp(a->resname, "PRP") )  {
      a->resnum = propane_num;
      
      if ( (propane_count%11) == 0 )
	propane_num++;
      
      propane_count++;
    }

    /* making sure resnum are sequential for isopropanol */
    else if ( !strcmp(a->resname, "POH") )  {
      a->resnum = isopropanol_num;
      
      if ( (isopropanol_count%12) == 0 )
	isopropanol_num++;
      
      isopropanol_count++;
    }
    /* making sure resnum are sequential for acetamide */
    else if ( !strcmp(a->resname, "ATM") )  {
      a->resnum = acetamide_num;
      
      if ( (acetamide_count%9) == 0 )
	acetamide_num++;
      
      acetamide_count++;
    }
    /* making sure resnum are sequential for isobutylene */
    else if ( !strcmp(a->resname, "IBE") )  {
      a->resnum = isobutylene_num;
      
      if ( (isobutylene_count%12) == 0 )
	isobutylene_num++;
      
      isobutylene_count++;
    }
    /* making sure resnum are sequential for benzene */
    else if ( !strcmp(a->resname, "BNZ") )  {
      a->resnum = bnz_num;
      
      if ( (bnz_count%12) == 0 )
	bnz_num++;
      
      bnz_count++;
    }
    /* making sure resnum are sequential for butane */
    else if ( !strcmp(a->resname, "BUT") )  {
      a->resnum = but_num;
      
      if ( (but_count%14) == 0 )
	but_num++;
      
      but_count++;
    }
    /* making sure resnum are sequential for GuH */
    else if ( !strcmp(a->resname, "GUH") )  {
      a->resnum = guh_num;
      
      if ( (guh_count%10) == 0 )
	guh_num++;
      
      guh_count++;
    }
    /* making sure resnum are sequential for methane */
    else if ( !strcmp(a->resname, "MTH") )  {
      a->resnum = mth_num;
#ifdef OFF      
      if (a->name[0] == 'C')   
	strncpy(a->name, "C4 ", 4);
#endif 
      if ( (mth_count%5) == 0 )
	mth_num++;
      
      mth_count++;
    }


    /* numbering the chlorines */
    else if ( STREQ(a->resname, "CL ") ) 
      a->resnum = cl_count_num++;


    /* numbering the sodiums */
    else if ( STREQ(a->resname, "NA ") ) 
      a->resnum = na_count_num++;


    /* the end of the protein */ 
    else if ( STREQ(a->resname, "END")) {
      strcpy(a->resname, f->atoms[atmcount-1]->resname);
      strcpy(a->name, "OXT");
      a->resnum = f->atoms[atmcount-1]->resnum;
    }
    
      
    /* deselecting all the hydrogens */
    a->restype = find_residue_type(a->resname) ; 
    a->definition = find_atom_definition(a->restype,a->name) ;
    
    if (a->name[0] == 'H') {
      SWITCH_ON(a,FLAG_1_DESELECT);
      SWITCH_ON(a,FLAG_2_DESELECT);
    } else {
      SWITCH_OFF(a,FLAG_1_DESELECT);
      SWITCH_OFF(a,FLAG_2_DESELECT);
    }
  }
}