#include "pdb.h"

/* incorporates modifications by Vossman on 11.17.2000
   to accomodate nucleic acids
*/

residue_description ** residue_types = NULL ;

int residue_type_number  ; 
int max_resuidue_type ; /* not realy maximum, because the program
				 will increment it as needed */

/* 
   find atom description inside a residue description. This
   canonicalize the name,so it can be a C string, i.e. null
   terminated.

   */

atom_definition * find_atom_definition_1(residue_description *residue, char *name)
{
	int atomnum = residue->atomnum ;
	atom_definition *ads =  &residue->atom_definitions[0];
	int index = 0 ;
        char first, second, third, fourth;
          first = name[0]; 
          second = name[1];
          third = name[2];
          fourth = name[3];
        if (first != '\0')
        {if (second == '\0')
         { second = ' ' ; third = ' '; fourth = '\0';}
         else if (third == '\0')
         { third = ' ' ; fourth = '\0' ; }
         else
           fourth='\0';
        }

        for(index= 0  ; index < atomnum ; index++)
	   {
	     atom_definition *ad = &ads[index] ;
		char *compare = ad->name;
                if (first == compare[0] && second == compare[1] &&
	            third == compare[2])
		  if(ad->flag)
		    return &ads[ad->flag] ;
		  else return ad ;
	   }
	return NULL ;
}
int handle_missing_atom_definition = 2 ;
residue_description *dumm_residue = NULL ;
atom_definition  *default_atom  = NULL;
atom_definition  *default_C_atom  = NULL;
atom_definition  *default_O_atom  = NULL;
atom_definition  *default_N_atom  = NULL;

atom_definition *assign_dummy_atom_definition(char * name)
{
  switch(name[0])
    { case 'N' : return default_N_atom ; break ;
       case 'C' : return default_C_atom ; break ;
       case 'O' : return default_O_atom ; break ;
    }
  return default_atom ;
}

atom_definition * find_atom_definition(residue_description *residue, char *name)
{
	atom_definition *ad ;
 	if (residue != NULL)
        {  
           ad = find_atom_definition_1(residue,name) ;
	   if (ad != NULL) return ad ;
	}
	if (dumm_residue != NULL) 
		ad = find_atom_definition_1 (dumm_residue, name) ;
	if (ad ==NULL)
	  switch (handle_missing_atom_definition)
	    {
	    case 0 :break ;
	    case 1 : if (residue == NULL ||residue->index >= 20) ad = assign_dummy_atom_definition(name) ;break ;
	    case 2 : ad = assign_dummy_atom_definition (name) ; break ;
	    }
	  if(ad == NULL)
	    serious_error(1, "Failed to find atom deinition",name);
	return ad ;
}
      

/* Find the type of the residue.  Return NULL if it is unknown. */

residue_description *find_residue_type (char *string)
{

  int ii ;
  char first, second, third;

  if(string[1] == '\0')
  {   first = ' ';
      second = ' ';
      third = string[0];
  } else if(string[2] == '\0')
  {   first = ' ';
      second = string[0];
      third = string[1];
  } else
  {   first = string[0];
      second = string[1];
      third = string[2];
  }
  for (ii= 0 ; ii< residue_type_number; ii++) {
    char *name = residue_types[ii]->name ;
    if( first == name[0] && second == name[1] && third == name[2] ) {
      residue_description *res =  residue_types[ii];
      return residue_types[res->index ] ; /* this is allowing redirection */
    }
  }
  
  return NULL; /* ??????? */
}



residue_description * alloc_residue (char *name, int atomnum, int dumm)
{


	if (residue_types == NULL)
	   {
		residue_types = (residue_description **)malloc_and_zero (5 * sizeof(char*));
		 max_resuidue_type = 5 ;
		 residue_type_number = 0 ;
	   }
     {
	int index ; 	
	residue_description *entry  = find_residue_type(name) ;

	if (NULL ==  entry)
	  {
		if (atomnum == 0 && !dumm) 
		  serious_error(1,"Trying to modify non-existing residue", name);

		
		if (residue_type_number >= max_resuidue_type) 
		   {
			int len =  max_resuidue_type * sizeof(char *) ;
			residue_description **new = (residue_description **)malloc_and_zero (len + len) ;
			
			bcopy (residue_types, new, len);
			free(residue_types) ;
			residue_types = new ;
			max_resuidue_type += max_resuidue_type;
		   }
		index = residue_type_number ;
		residue_type_number++ ;
	   }
	else
	  {
		index = entry->index ;
		if (atomnum != 0 && atomnum != entry->atomnum)
		    {
			free(entry) ;
			entry = NULL ;
		   }
	   }
	if (entry == NULL)
	     {
		int entry_len = sizeof(residue_description) + (atomnum -1) * sizeof(atom_definition) ; 
		entry = (residue_description *) malloc_and_zero(entry_len);
		residue_types[	index ] = entry ;
		strncpy(entry->name, name,3) ;
		entry->atomnum = atomnum ;
		entry->index = index  ;
             }
	return entry ;
      }
}
	


/*  read a residue description. The number of atoms, name and index is already set. */

void read_residue_atoms (FILE *ff, residue_description * entry)
{
	int ii ,type;
	char line[200] ;
	char string [200] ; 
	char type_string [200] ; 
	double vdw, covalent , charge;


	while (fgets(line,200,ff) != NULL)	
	    {
/* ATOM NAME */
	 	ii = sscanf(line,LF_STRING_06, string, &type_string, &vdw,&covalent,&charge) ;
		if (ii <= 0 ) break ;

                if(string[1] == '\0' || string[1] == ' ')
                          { string[1] = ' ' ; string[2] = ' ' ;string[3] = '\0';}
                else if(string[2] == '\0' || string[1] == ' ')
                          { string[2] = ' ' ; string[3] = '\0';}
                else if(string[3] == '\0') ;
                else if(string[4] == '\0') ;
                else serious_error(1, "Wrong Atom name (too long)", string) ;

                /* NO SHIFTING
		if(string[1] != '\0' && string[2] != '\0' &&
                        string[3] != '\0' && string[4] != '\0')
		    serious_error(1, "Wrong Atom name (too long)", string) ;*/
	    {
		atom_definition *ad =  find_atom_definition_1(entry, string) ;
		if (ad == NULL)
		  {  	char search[4] = {'\0','\0','\0','\0'} ;
			ad = find_atom_definition_1(entry,search) ;
			if (ad == NULL)
				serious_error(1,"Too many atoms in residue definition", entry->name);
			strncpy(ad->name, string, 4) ;
		  }
		if (ii > 1)
                   if (type_string[0] == '=')
		     { 	
                        /*DO NOT need new shifted string*/
                        { 
                          atom_definition * real_atom = find_atom_definition(entry,type_string+1);
			  if (real_atom == NULL)
			    serious_error(1, "Failed to find real atom for dummy atom", string) ;
			  ad->vdw = real_atom->vdw;
			  ad->covalent = real_atom->covalent;
			  ad->covalent = real_atom->covalent;
			  ad->flag =  PROTEIN_ATOM_DEFINITION_INDEX(real_atom,entry) ;
                        }
                      }
                  else
		    {
		      	int type = find_type_definition (type_string) ; 
		  
			ad->type_index = type ;
			ad->flag = 0 ;
		   	if (ii > 2 && vdw != 0.0)
		  	    ad->vdw = vdw ;
		   	else if (type != 0)
			    ad->vdw = atom_vdw[type] ;

		   	if (ii > 3 && covalent != 0.0)
		  	   ad->covalent = covalent ;
		   	else if (type != 0)
			   ad->covalent = atom_covalent[type] ;

		   	if (ii > 4 && charge != 0.0)
		  	    ad->charge = charge ;
		   	else if (type != 0)
			       ad->charge = atom_charge[type] ;
		  }
		
	   }}
}
			
		  


/* read description of residues */

void read_residue_definitions (FILE *ff)
{
	int ii,atomnum,kk;
	char line[200] ; 
	char string [200] ; 

	while (fgets(line,200,ff) != NULL)
	   {
/* RESI NAME */
	 	ii = sscanf(line,"%s%d", string, &atomnum) ;

		if(ii > 0 && string[0] != '!' && string[0] != ';') 
		    {	
			int real_index = -1 ; 

                /* A fix for simalarity to pdb files*/
                        if(string[1] == '\0')
                            { string[1] = ' ' ; string[2] = string[0] ;
                              string[0] = ' ' ; }
                        else if(string[2] == '\0')
                            { string[2] = string[1] ; string[1] = string[0] ;
                              string[0] = ' ' ; }
                        else if(string[3] != '\0')
			  serious_error(1, "Wrong residue name (too long)", string) ;
                        string[3] = '\0';

                        if  (ii == 1)
			  {
                            char *poi = line ;
			    
			    while(isspace(*poi)) poi ++;
			    while(!isspace(*poi)) poi ++;
			    while(isspace(*poi)) poi ++;
                            
/*Alias does not work for nucleic acids!!!*/
			    if (*poi == '=')
			      {
				residue_description *alias = find_residue_type(poi+1) ;

				if (alias == NULL)
				  serious_error(ARG_TYPE_STRING, "Failed to find alias residue", line) ;
				real_index = alias->index ;
				atomnum = 0 ;
			      }
			    else atomnum = - 1 ;
			  }

			if ( 0 > atomnum || atomnum > 200) 
			   serious_error(1, "Wrong atom number in residue type", string) ;
		
			{
		  	   residue_description *entry = alloc_residue(string, atomnum , real_index >= 0 ? 1 : 0) ;
			   if (real_index >= 0)
			     entry->index = real_index ; 
		  	   read_residue_atoms (ff,entry) ;
			}
		    }
	   }
}
		    












/*   traverse all the atom definitions in a residue  */
void traverse_residue_definition (residue_description *res, void (*function)(), void *arg)
{
    int atomnum = res->atomnum ;
    atom_definition * ads = res->atom_definitions ;
    int ii ;

   for (ii = 0 ; ii < atomnum ; ii++)
     {
	if(ads->name[0] != '\0' )
		(*function)(ads, arg) ;
	ads++ ;
    }
}

void traverse_all_residue_definitions (void (*function)(), void *arg, void *arg1)
{
   int ii ;
  
   for(ii = 0 ; ii< residue_type_number ; ii++)
         (*function)(residue_types[ ii ] , arg,arg1) ;
}


void traverse_atom_definitions (void (*function)(), void *arg)
{
   traverse_all_residue_definitions(  traverse_residue_definition, function,arg) ;
  /* traverse_additional_atom_definitions */
}

void read_standard_volumes(char *filename)
{
  FILE *ff = fopen(filename,"r") ;
  char string[200] ;
  int ii ;
  residue_description *res ; 
  int ssbond = 0 ;
  for(ii = 0 ; ii < 21 ; ii++)
    {
      while(1)
	{
	  fgets(string,200, ff) ;
	  if (isupper(string[0]) && isupper(string[1]) && isupper(string[2]))
	    break;
	}
      res = find_residue_type(string) ;
      while(1)
	{
	  int count;
	  char name[10];
	  double volume ;
	  fgets(string, 200, ff) ;
/* ATOM NAME */
	  if (sscanf(string,"%s %d %lf", name, &count,&volume) < 3)  break ;
	  {
	    atom_definition *atom = find_atom_definition(res, name) ;
	    if (ssbond)
	      atom->eradsq = volume;
	    else atom->eradius = volume ;
	  }
	}
      if(!strcmp (res->name, "CYS") && !ssbond)
	ssbond  = 1;
      else ssbond = 0 ;
    }
}

/*
New routine to handle the reading of both protein and nucleic acid 
standard volumes 
*/

void read_standard_volumes_nt(char *filename)
{
  FILE *ff = fopen(filename,"r") ;
  char string[200] ;
  int ii ;
  residue_description *res ; 
  int ssbond = 0 ;
  for(ii = 0 ; ii < 25 ; ii++) /* 25 = 21 aa + 4 na */
    {
      do
	{
	  fgets(string,200, ff) ;
	} while( !isupper(string[0]) );
/* RESI NAME */
      string[3] = '\0';
      if(string[1] == '\0' || string[1] == '\n' || string[1] == ' ')
      { string[1] = ' ' ; string[2] = string[0] ;
        string[0] = ' ' ; string[3] = '\0'; }
      else if(string[2] == '\0' || string[2] == '\n' || string[2] == ' ')
      { string[2] = string[1] ; string[1] = string[0] ;
        string[0] = ' ' ; string[3] = '\0'; }
      else if(string[3] != '\0'  || string[3] == '\n' || string[3] == ' ')
	serious_error(1, "Wrong residue name (too long)", string) ;
      string[3] = '\0';
      res = find_residue_type(string) ;
      while(1)
	{
	  int count;
	  char name[10];
	  double volume ;
	  fgets(string, 200, ff) ;
/* ATOM NAME */
	  if (sscanf(string,"%s %d %lf", name, &count,&volume) < 3)  break ;
	  {
	    atom_definition *atom = find_atom_definition(res, name) ;
	    if (ssbond)
	      atom->eradsq = volume;
	    else atom->eradius = volume ;
	  }
	}
      if(!strcmp (res->name, "CYS") && !ssbond)
	ssbond  = 1;  
      else
	ssbond = 0 ;
    }
}


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


/*   ===================                                         */
void ShowAllResidueTypes  
/*   ===================                                         */
/* (int residue_type_number, residue_description ** residue_types ) */
(void)
{
  int ii ;
  for (ii= 0 ; ii< residue_type_number; ii++) {
    residue_description * r = residue_types[ii] ;
    printf("@ %8x  %3d : %3s = %4d (%4d atoms)\n",
	   residue_types[ii],ii,r->name,r->index,r->atomnum);
  }
}