/*                                                  +++++++++          */
# define FILENAME                                   "pdbio.c"          
/*                                                  +++++++++          */
# include "util.h"

/*  =======                                                       */
int nrd_pdb
/*  =======                                                       */
  (FILE *stream, PDB_data *P, double xx[], double yy[], double zz[])
{
  P->natom = rd_pdb(stream, P->iser, P->atnam,P->resnam,P->chnnam,P->iresn,
		xx,yy,zz,P->occ,P->b,P->segid);
  return P->natom;
} 

/*  ======                                                       */
int rd_pdb
/*  ======                                                       */
	  (FILE *stream, int iser[], char atnam[][5], char resnam[][5], 
	   char chnnam[], int iresn[], double xx[], double yy[], double zz[], 
	   double occ[], double b[], char segid[][5])
{
  int i;
  static char buf[BUFSIZ+1];
  char *buf1 = buf + 1;
  int natom = 0;

  if (stream == NULL) 
    nonfatal(__FILE__,"rd_pdb","stream == NULL");
  while ((fgets(buf1,BUFSIZ,stream) != NULL) &&
	 (!STRNEQ(buf1,"END",3))) {
#   ifdef OFF
    printf("DEBUG: buf1=%s\n",buf1);
#   endif
    if ( STRNEQ(buf1,"ATOM",4) || STRNEQ(buf1,"HETATM",6) ) { 
#     ifdef OFF
      printf("DEBUG: buf1=%s\n",buf1);
#     endif
      if (++natom > MXPDB) {
	printf ("DEBUG: natom=%4d MXPDB=%4d\n",natom,MXPDB);
	nonfatal(FILENAME,"rd_pdb","natom > MXPDB");
      }
      sscanf(buf+7,"%5d",iser+natom);
      atnam[natom][0]=buf[13];
      atnam[natom][1]=buf[14];
      atnam[natom][2]=buf[15];
      atnam[natom][3]=buf[16];
      atnam[natom][4]='\0';
      resnam[natom][0]=buf[18];
      resnam[natom][1]=buf[19];
      resnam[natom][2]=buf[20];
      resnam[natom][3]=buf[21];
      resnam[natom][4]='\0';
      chnnam[natom]=buf[22];
      sscanf(buf+23,"%4d",iresn+natom); 
      sscanf(buf+31,"%8lf",xx+natom);
      sscanf(buf+39,"%8lf",yy+natom);
      sscanf(buf+47,"%8lf",zz+natom);
      sscanf(buf+55,"%6lf",occ+natom);
      sscanf(buf+61,"%6lf",b+natom);
      segid[natom][0]=buf[73];
      segid[natom][1]=buf[74];
      segid[natom][2]=buf[75];
      segid[natom][3]=buf[76];
      segid[natom][4]='\0';
#     ifdef OFF
      pdbdebug(iser, atnam, resnam,chnnam, iresn, xx, yy, 
	       zz, occ, b,segid,natom,buf1);
#     endif
    }
  }
  return(natom);
}

/*   ========                                                           */
void pdbdebug
/*   ========                                                           */
	  (int iser[], char atnam[][5], char resnam[][5], 
	   char chnnam[], int iresn[], double xx[], double yy[], double zz[], 
	   double occ[], double b[], char segid[][5], int natom, char *buf)
{
  int i=natom;
  printf("buf: %s\n",buf);
  printf(PDBFMT,iser[i],atnam[i],resnam[i],chnnam[i],iresn[i],
	  xx[i],yy[i],zz[i],occ[i],b[i],segid[i]);
}  

/*   ===========                                                       */
void nprt_pdbrec
/*   ===========                                                       */
  (PDB_data *P,double xx[], double yy[], double zz[],int i)
{
  printf(PDBFMT,P->iser[i],P->atnam[i],P->resnam[i],P->chnnam[i],P->iresn[i],
	  xx[i],yy[i],zz[i],P->occ[i],P->b[i],P->segid[i]);
}  

/*   ========                                               */
void nwrt_pdb
/*   ========                                               */
	  (FILE *stream, PDB_data *P, 
	   double xx[], double yy[], double zz[])
{
  int i = 0;
  while (++i<=P->natom) 
    fprintf(stream,PDBFMT,P->iser[i],P->atnam[i],P->resnam[i],P->chnnam[i],
	    P->iresn[i],xx[i],yy[i],zz[i],P->occ[i],P->b[i],
	    P->segid[i]);
  
}

/*   ========                                               */
void nwrt_dia
/*   ========                                               */
	  (FILE *stream, PDB_data *P, 
	   double xx[], double yy[], double zz[])
{
  char atm, res[5];
  int k, i = 0;

  fprintf(stream,
	  "    1.000     1.000     1.000    90.000    90.000    90.000\n");
  fprintf(stream,
	  "\n(4F10.5, 3I5, F9.2, 1X, 2A4, 3X, A4)\n");


/* 
** Determine the atom type number
*/
    while (++i<=P->natom) {
      atm = P->atnam[i][1];
      if (atm == 'C')
	k = 1;
      else if (atm == 'O')
	k = 3;
      else if (atm == 'N')
        k = 2;
      else if (atm == 'S')
        k = 4;
      else if (atm == 'P')
        k = 5;
      else if (atm == 'M')
	k = 6;

      sprintf(res,"%d",P->iresn[i]);
      fprintf (stream,DIAFMT,
	       xx[i],yy[i],zz[i],P->b[i],
	       k,P->iresn[i],i,
	    P->occ[i],
	       P->resnam[i],res,P->atnam[i]);
    }
}

/****
/ATOM|HETATM/ {
 Is there a chain identifier?
   chain = 0
   if (substr($0,22,1) != " ") {
      chain = 1
      str = $5$6
   }
.
.
.
 Are we starting a new residue?
   if (str != oldres) {
      res++
      oldres = str
   }
 Print this residue in Diamond format.
   if (chain == 1) {
      printf (diafmt,$7,$8,$9,$11,k,res,++j,$10,$4,str,$3)
          x  y  z  b  ityp ires iatom restyp resnam atmnam
   }
   else {
****/


/*   =======                                                   */
void wrt_pdb
/*   =======                                                   */
	  (FILE *stream, int iser[], char atnam[][5], char resnam[][5], 
	   char chnnam[], int iresn[], double xx[], double yy[], double zz[], 
	   double occ[], double b[], char segid[][5], int natom)
{
  int i;
  for(i=1;i<=natom;i++) {
    fprintf(stream,PDBFMT,iser[i],atnam[i],resnam[i],chnnam[i],iresn[i],
	    xx[i],yy[i],zz[i],occ[i],b[i],segid[i]);
  }
}

/*   ===========                                               */
void nwrt_pdbsel
/*   ===========                                               */
	  (FILE *stream, PDB_data *P, 
	   double xx[], double yy[], double zz[], int sel[])
{
  int i = 0;
  while (++i<=P->natom) 
    if (sel[i]) 
      fprintf(stream,PDBFMT,P->iser[i],P->atnam[i],P->resnam[i],P->chnnam[i],
	      P->iresn[i],xx[i],yy[i],zz[i],P->occ[i],P->b[i],
	      P->segid[i]);

}

/*   ==========                                               */
void wrt_pdbsel
/*   ==========                                               */
	  (FILE *stream, int iser[], char atnam[][5], char resnam[][5], 
	   char chnnam[], int iresn[], double xx[], double yy[], double zz[], 
	   double occ[], double b[], char segid[][5], int natom, int sel[])
{
  int i = 0;
  while (++i<=natom) 
    if (sel[i]) 
      fprintf(stream,PDBFMT,iser[i],atnam[i],resnam[i],chnnam[i],iresn[i],
	      xx[i],yy[i],zz[i],occ[i],b[i],segid[i]);
}

/*   =============                                               */
int  SimpleWrtAtom
/*   =============                                               */
  (FILE *stream, int iresn, double v[3])
{
  static int 
    iser = 0;
  char * segid = "     " ;
  char chnnam = ' ' ;
  char * atnam = " CA " ;
  char * resnam = "UNK";
  double occ = 1 ;
  double b = 0;
  iser++;
  fprintf(stream,PDBFMT,iser,atnam,resnam,chnnam,iresn,
	  v[1],v[2],v[3],occ,b,segid);
  return iser ;
}

/*   ==============                                               */
int  SimpleWrtAtom2
/*   ==============                                               */
  (FILE *stream, int iresn, char * resnam, double v[3])
{
  static int 
    iser = 0;
  char * segid = "     " ;
  char chnnam = ' ' ;
  char * atnam = " CA " ;
  double occ = 1 ;
  double b = 0;
  iser++;
  fprintf(stream,PDBFMT,iser,atnam,resnam,chnnam,iresn,
	  v[1],v[2],v[3],occ,b,segid);
  return iser ;
}



/*  =======                                                               */
int Nselect (PDB_data *P,int sel[])
/*  =======                                                               */
{
  char es; 
  int i, nsel = 0;
  
  es = *getenv("PDB_SELECT");
  for(i=1;i<=P->natom;i++) 
    if 
      (es == 'O' 
	&& P->atnam[i][1]=='O' && P->atnam[i][2]=='H' && P->atnam[i][3]=='2' ) {
      sel[i]=1;
      nsel++;
    } else if 
      (es == 'H' 
       && P->atnam[i][1]=='H' 
       && (P->atnam[i][2]=='1' || P->atnam[i][2]=='2' )) {
      sel[i]=1;
      nsel++;
    } else {
      sel[i]=0;
    }
  printf("%s (select): es=%c gives nsel=%d from P->natom=%d\n",
	 FILENAME,es,nsel,P->natom);
  return nsel;
}
 
/*
** Positional selection:
** 
    for(i=1;i<=dcd->natom;i++) 
      if (fabs(xx[i]) < 8.0 && fabs(yy[i]) < 8.0 && fabs(zz[i])<8.0)
	sel[i] = 1;
      else
	sel[i] = 0;
**
*/	  

/*       -------                                                        */
# define ISBBONE(atnam,resnam)\
/*       -------                                                        */\
 ((STREQ(atnam," CA ") || STREQ(atnam," O  ") ||\
   STREQ(atnam," C  ") || STREQ(atnam," N  ")) &&\
  !(STREQ(resnam,"WAT ") || STREQ(resnam,"TIP3")))
			 
/*   ==========                                                  */
void pdb_select(int argc, char *argv[])
/*   ==========                                                  */
{
  bool 
    bbonly = (argc >= 3  ? STREQ(argv[3],"-bb") : 0 );
  PDB_data pdb, *P = &pdb;
  double xx[MXPDB],yy[MXPDB], zz[MXPDB];
  int 
    sel[MXPDB],i,
    ires_min = atoi (argv[1]),
    ires_max = atoi (argv[2]),
    natom = nrd_pdb(stdin,P,xx,yy,zz);
  for (i=1; i<natom; i++) {
    sel[i] = (P->iresn[i] >= ires_min && P->iresn[i] <= ires_max );
    if (bbonly && !ISBBONE(P->atnam[i],P->resnam[i]))
      sel[i] = 0;
  }
  nwrt_pdbsel(stdout,P,xx,yy,zz,sel);
} 
# undef ISSBONE

/*   ================                                                     */
void pdb2vectcentroid (int argc, char *argv[])
/*   ================                                                     */
{
  double xx[MXPDB],yy[MXPDB], zz[MXPDB];
  PDB_data pdb, *P = &pdb;
  int 
    i, 
    movedraw = (argc >= 1 ? atoi(argv[1]) : 1),
    natom = nrd_pdb(stdin,P,xx,yy,zz);
  Vector v;

  HERE_OFF("pdb2vectcentroid BEGINS");
  for (i=1; i<=natom; i++) {
    v[0] += xx[i];
    v[1] += yy[i];
    v[2] += zz[i];
  }
  v[0] /= natom; v[1] /= natom; v[2] /= natom;
  WRT1VECT(stdout,v,movedraw);
}