#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <string.h>
#include <stddef.h>
#include <ctype.h>
#include "io_lib_header.h"
#include "util_lib_header.h"
#include "dp_lib_header.h" 
#include "define_header.h"



#define FATAL "fatal:reformat_struc"

char * process_repeat (char *aln, char *seq, char *pdb)
{
  char *tf, *file, *name;
  Alignment *A, *A2;
  Sequence *S, *P;
  int r1, r2, is_r1, is_r2, l1=0, l2=0, a;
  int *pos;
  FILE *fp;
  
  A=main_read_aln (aln, NULL);
  for (a=0; a<A->nseq; a++)ungap(A->seq_al[a]);
  
  S=main_read_seq (seq);
  P=get_pdb_sequence (pdb);
  

  A2=align_two_sequences (S->seq[0], P->seq[0], "pam250mt", -10, -1, "myers_miller_pair_wise");
  
  pos=vcalloc ( A2->len_aln+1, sizeof (int));

  for (l1=0, l2=0,a=0; a< A2->len_aln; a++)
    {
      
      r1=A2->seq_al[0][a];
      r2=A2->seq_al[1][a];
      
      is_r1=1-is_gap(r1);
      is_r2=1-is_gap(r2);
      
      l1+=is_r1;
      l2+=is_r2;
      
      if (!is_r1);
      else
	{
	  pos[l1]=l2;
	}
    }
  tf=vtmpnam(NULL);
  fp=vfopen (tf, "w");
  for (a=0; a<A->nseq; a++)
    {
      int *coor, b, c;
      
      name=A->name[a];
      file=vtmpnam (NULL);
      coor=string2num_list2 (name, "-");

      //Check the compatibility between the guide sequence and the coordinates
      for ( c=0,b=coor[1]-1; b<coor[2]; b++, c++)
	{

	  if (tolower(A->seq_al[a][c])!=tolower(S->seq[0][b]))
	    printf_exit (EXIT_FAILURE, stderr, "Incompatibility between the repeat [%s] and the master Sequence [%s]\n%s",A->name[a], seq, A->seq_al[a]);
	}
   
      printf_system ( "extract_from_pdb %s -coor %d %d -seq_field SEQRES> %s", pdb,pos[coor[1]-1],pos[coor[2]-1], file);
      fprintf (fp, ">%s _P_ %s\n", name, file);
      vfree (coor);
    }
  vfclose (fp);
  return tf;
}
  
  
  
char *     normalize_pdb_file  (char *name_in, char *seq, char *out_file)
{
  char command[1000];
  Sequence *S;
  Alignment *A;
  int a;
  int start, end, npdb, r1, r2;
  char name[100];

  
 if ( !name_in) return NULL;
 else
   {
     sprintf ( name, "%s", name_in);
   }
 
 if ( !is_pdb_file(name))
    {
      fprintf(stdout, "\nERROR[normalize_pdb_file]: %s is not a pdb file[FATAL:%s]\n", name, PROGRAM);
      myexit (EXIT_FAILURE);
    }

  S=get_pdb_sequence (name);
  A=align_two_sequences (S->seq[0],seq,"idmat",-3,0, "fasta_pair_wise");
 
  

  for (start=-1, end=-1,npdb=0,a=0; a< A->len_aln; a++)
    {
      r1=1-is_gap(A->seq_al[0][a]);
      r2=1-is_gap(A->seq_al[1][a]);
      
      npdb+=r1;

      if (r1 && r2 && start==-1)start=npdb;
      if (r1 && r2)end=npdb;
    }

  free_aln(A);  
  free_sequence (S, -1);
  
  sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -coor %d %d -nodiagnostic > %s", check_file_exists(name), start, end, out_file);
  my_system ( command);
  return out_file;
  }

Ca_trace * trim_ca_trace (Ca_trace *T, char *seq )
{
  /*This function:
    -removes from Ca trace all the residues that are not in the sequence
    -add in the Ca trace residues unmatched in the structure (it gives them a NULL structure)
  */
  Alignment *ALN;
  Atom *A;  
  int a,l, s, r, is_r, is_s;
  int *seq_cache, *struc_cache;
  int tot_l=0;

  char buf1[10000];
  char buf2[10000];

  /*  lower_string (T->seq);
      lower_string (seq);
  */


  sprintf (buf1, "%s", T->seq);
  sprintf (buf2, "%s", seq);
  lower_string (buf1);
  lower_string (buf2);
  

  if ( strm (buf1,buf2))return T;
  else
    {
     ALN=align_two_sequences (T->seq,seq, "est_idmat",-1, 0,"fasta_pair_wise"); 
     struc_cache=vcalloc (ALN->len_aln+1, sizeof (int));
     seq_cache  =vcalloc (ALN->len_aln+1, sizeof (int));
     
     for ( r=0, s=0,a=0; a< ALN->len_aln; a++)
       {
	 is_r=!is_gap(ALN->seq_al[0][a]);
	 is_s=!is_gap(ALN->seq_al[1][a]);
	 
	 r+=is_r;
	 s+=is_s;
	 
	 if ( is_s && is_r)
	   {
	     struc_cache[r-1]=s-1;
	     seq_cache[s-1]=r-1;
	   }
	 else if ( is_s && !is_r)
	   {
	     seq_cache[s-1]=-1;
	     
	   }
	 else if ( !is_s && is_r)
	   {
	     struc_cache[r-1]=-1;
	   }
       }
     
      T->ca=vrealloc ( T->ca, sizeof (Atom*)*(ALN->len_aln+1));
      T->peptide_chain=vrealloc ( T->peptide_chain, (sizeof (Amino_acid*))*(ALN->len_aln+1));
      T->seq=vrealloc ( T->seq,   ALN->len_aln+1);

      for ( a=T->len; a< ALN->len_aln; a++)
	{
	  T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));
	}
      
      
      /*Read every atom*/
      for ( a=0; a< T->n_atom; a++)
	{
	  
	  A=(T->structure[a]);
	  if ( struc_cache[A->res_num-1]==-1)continue;
	  else
	    {
	      /*set the struc residue to its sequence index*/
	      A->res_num=struc_cache[A->res_num-1]+1;
	      if (strm (A->type, "CA")) {T->ca[A->res_num-1]=A;tot_l++;}
	      if ( strm (A->type, "CA"))(T->peptide_chain[A->res_num-1])->CA=A;
	      if ( strm (A->type, "C"))(T->peptide_chain[A->res_num-1] )->C=A;
	      if ( strm (A->type, "CB"))(T->peptide_chain[A->res_num-1])->CB=A;
	      if ( strm (A->type, "N"))(T->peptide_chain[A->res_num-1] )->N=A;
	    }
	}
 
      l=strlen(seq);     
      for ( a=0;a< l; a++)
       {

	 if ( seq_cache[a]==-1)
	   {
	     tot_l++;
	     T->ca[a]=NULL;
	     
	     if (!T->peptide_chain[a])T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));
	     T->peptide_chain[a]->CA=NULL;
	     T->peptide_chain[a]->C =NULL;
	     T->peptide_chain[a]->CB=NULL;
	     T->peptide_chain[a]->N=NULL;
	     T->seq[a]='x';
	   }
	 else
	   {
	     T->seq[a]=seq[a];
	   }
       }
      T->len=ALN->len_aln;
     
      /*
	T->len=tot_l;
      */
      
     free_aln (ALN);
     vfree(seq_cache);
     vfree(struc_cache);
     
	  
    }
  return T;
}

Ca_trace * read_ca_trace (char *name, char *seq_field )
{
  char *tp_name=NULL;
  char command[10000];
  

  if ( !is_simple_pdb_file (name))
    {
      tp_name=vtmpnam (NULL);
      sprintf ( command, "extract_from_pdb -seq_field %s  -infile %s -atom ALL -chain FIRST -mode simple> %s",seq_field, check_file_exists(name), tp_name);
      if ( getenv4debug ("DEBUG_EXTRACT_FROM_PDB"))fprintf ( stderr, "\n[DEBUG_EXTRACT_FROM_PDB:read_ca_trace] %s\n", command);
      my_system ( command);
    }
  else
    tp_name=name;
  
  return simple_read_ca_trace (tp_name);
}

Ca_trace * simple_read_ca_trace (char *tp_name )
    {
	/*This function reads a pdb file into a Ca_trace structure*/
	 
	 

	int a, c, n;
	FILE *fp;
	Atom *A;
	char res;
	char *buf;
	Ca_trace *T=NULL;
	int res_num=0, last_res_num=0;

	
	buf=vcalloc ( VERY_LONG_STRING, sizeof (char));
	n=count_n_line_in_file (tp_name );
	
	if ( !T)
	    {
	    T=vcalloc ( 1, sizeof ( Ca_trace));
	    declare_name (T->name);
	    }
	
	/*1 Get the complete sequence: replace missing residues with Xs*/
	for (a=0; a< VERY_LONG_STRING; a++)buf[a]='x';
	res=res_num=0;
	
	fp=vfopen (tp_name, "r");
        while ( (c=fgetc(fp))!='>');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');	
	while ((c=fgetc(fp))!=EOF)
	  {
	       ungetc(c, fp);

	       fscanf (fp, "%*s %*s %c %*c %d %*f %*f %*f\n",&res,&res_num);
	       if (res)
		 {
		  
		   res=tolower (res);
		   buf[res_num-1]=res;
		   last_res_num=res_num;
		 }
	       res=res_num=0;
	  }
	buf[last_res_num]='\0';
	vfclose (fp);
	/*Sequence Read*/
	

	T->len=strlen (buf);
	T->seq=vcalloc ( T->len+1, sizeof (char));
	buf=lower_string (buf);
	sprintf ( T->seq, "%s", buf);
	n+=T->len;
	T->structure=vcalloc ( n, sizeof (Atom*));
	for ( a=0; a< n; a++)T->structure[a]=vcalloc ( 1, sizeof (Atom));
	T->ca=vcalloc ( T->len+1, sizeof ( Atom*));
	a=0;

	fp=vfopen (tp_name, "r");
        while ( (c=fgetc(fp))!='>');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');
        fscanf ( fp, "%*s" );
        while ( (c=fgetc(fp))!='\n');	
	
	while ((c=fgetc(fp))!=EOF)
	   {
	       ungetc(c, fp);
	       A=T->structure[a];
	       A->num=a;
	       fscanf (fp, "%*s %s %s %*c %d %f %f %f\n",A->type, A->res,&A->res_num, &A->x, &A->y, &A->z);
	       res=A->res[0];
	       
	       res=tolower (res);
	       
	       T->seq[A->res_num-1]=res;

	       if ( strm ( A->type, "CA")) T->ca[A->res_num-1]=A;
	       a++;
	   }
	T->n_atom=a;

	T->peptide_chain=vcalloc (T->len+1, sizeof (Amino_acid*));
	for ( a=0; a<=T->len; a++) T->peptide_chain[a]=vcalloc (1, sizeof (Amino_acid));				   
	for ( a=0; a< T->n_atom; a++)
	  {
	    A=T->structure[a];

	    if ( strm (A->type, "CA"))(T->peptide_chain[A->res_num-1])->CA=A;
	    if ( strm (A->type, "C"))(T->peptide_chain[A->res_num-1] )->C=A;
	    if ( strm (A->type, "CB"))(T->peptide_chain[A->res_num-1])->CB=A;
	    if ( strm (A->type, "N"))(T->peptide_chain[A->res_num-1] )->N=A;
	  }
	    

	vfclose (fp);	
	vfree (buf);
	return T;
    }
Ca_trace * hasch_ca_trace    ( Ca_trace *T)
    {
      
      T=hasch_ca_trace_nb (T);
      T=hasch_ca_trace_bubble (T);
      T=hasch_ca_trace_transversal (T);
      return T;
    }
Ca_trace * hasch_ca_trace_transversal ( Ca_trace *TRACE)
    {
	/*This function gets the Coordinates of a protein and computes the distance of each Ca to its 
	  
	  given a Ca,
	  Compute the distance between, CA-x and CA+x with x=[1-N_ca]
	  T->nb[a][0]-->Number of distances.
	  T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
	  T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
	*/
	
	int a, b, d;
	float dist;
	Atom *A, *B;
		
	Struct_nb *T;
	Pdb_param *PP;


	TRACE->Transversal=vcalloc ( 1, sizeof (Struct_nb));
	
	T=TRACE->Transversal;
	PP=TRACE->pdb_param;

	if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
	if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);

	for (d=0,a=0; a< TRACE->len; a++)
	    {

	        for ( b=1; b<=PP->N_ca; b++)
		    {
		    if ( (a-b)<0 || (a+b)>=TRACE->len)continue;
		    A=TRACE->ca[a-b];
		    B=TRACE->ca[a+b];
		    dist=get_atomic_distance ( A, B);
		    
		    T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
		    T->nb[a][T->nb[a][0]]=b;

		    T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
		    T->d_nb[a][T->nb[a][0]]=dist;
		    
		    d++;
		    }
		T->max_nb=MAX (T->max_nb, T->nb[a][0]);
	    }
	return TRACE;

    }			       

Ca_trace * hasch_ca_trace_nb ( Ca_trace *TRACE)
    {
	/*This function gets the Coordinates of a protein and computes the distance of each Ca to its 
	  T->N_ca Ngb.
	  The first Ngb to the left and to the right are excluded
	  Ngd to the left get negative distances
	  Ngb to the right receive positive distances
	  T->nb[a][0]-->Number of ngb.
	  T->nb[a][1... T->nb[a][0]]-->ngb index with respect to the Ca chain
	  T->d_nb[a][1... T->d_nb[a][0]]-->ngb index with respect to the Ca chain
	*/
	


	int a, b, d;
	float dist;
	Atom *A, *B;
		
	Struct_nb *T;
	Pdb_param *PP;
	
	TRACE->Chain=vcalloc ( 1, sizeof (Struct_nb));
	
	T=TRACE->Chain;
	PP=TRACE->pdb_param;

	if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
	if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);

	for (d=0,a=0; a< TRACE->len; a++)
	    {
		for ( b=MAX(0,a-PP->N_ca); b< MIN( a+PP->N_ca, TRACE->len); b++)
		        {
			if (FABS(a-b)<2)continue;
			A=TRACE->ca[a];
			B=TRACE->ca[b];
			if ( !A || !B)continue;
			dist=get_atomic_distance ( A, B);
			if (b<a) dist=-dist;
		

			T->nb[a]=vrealloc ( T->nb[a], (++T->nb[a][0]+1)*sizeof (int));
			T->nb[a][T->nb[a][0]]=b;
			
			T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
			T->d_nb[a][T->nb[a][0]]=dist;
			d++;
			}
		    
		T->max_nb=MAX (T->max_nb, T->nb[a][0]);
	    }
	return TRACE;

    }	
Ca_trace * hasch_ca_trace_bubble ( Ca_trace *TRACE)
    {
      
	int a, b;
	float dist;
	Atom *A, *B;
	float **list;
	Pdb_param *PP;
	Struct_nb *T;

	PP=TRACE->pdb_param;
	TRACE->Bubble=vcalloc ( 1, sizeof (Struct_nb));
	T=TRACE->Bubble;
	
		
	
	if ( !T->nb)T->nb=declare_int (TRACE->len+1, 1);
	if ( !T->d_nb)T->d_nb=declare_float (TRACE->len+1, 1);
	list=declare_float ( TRACE->n_atom, 3);
		

	for (a=0; a< TRACE->len; a++)
	    {
		for ( b=0; b< TRACE->len; b++)
		    {
			A=TRACE->ca[a];
			B=TRACE->ca[b];
			if ( !A || !B)continue;
			dist=get_atomic_distance ( A, B);
			
			if ( dist<PP->maximum_distance && FABS((A->res_num-B->res_num))>2)
			   {
			     T->nb[a][0]++;			       			       
			     T->nb[a]=vrealloc ( T->nb[a], (T->nb[a][0]+1)*sizeof (int));
			     T->nb[a][T->nb[a][0]]=(TRACE->ca[b])->num;
			     
			     T->d_nb[a]=vrealloc ( T->d_nb[a], (T->nb[a][0]+1)*sizeof (float));
			     T->d_nb[a][T->nb[a][0]]= ((a<b)?dist:-dist);
			   }						
		    }
		T->max_nb=MAX (T->max_nb, T->nb[a][0]);
	
	    }
	
	for ( a=0; a< TRACE->len; a++)
	    {
	    for ( b=0; b< T->nb[a][0]; b++)
	        {
		    list[b][0]=T->nb[a][b+1];
		    list[b][1]=T->d_nb[a][b+1];
		    list[b][2]=(TRACE->structure[T->nb[a][b+1]])->res_num;
		}
	    
	    sort_float ( list, 3,2, 0, T->nb[a][0]-1);
	    for ( b=0; b< T->nb[a][0]; b++)
	        {
		    T->nb[a][b+1]=list[b][0];
		    T->d_nb[a][b+1]=list[b][1];
		}
	    }

	free_float ( list, -1);
	return TRACE;

    }

       
float ** measure_ca_distances(Ca_trace *T)
   {
       int a, b;
       Atom *A, *B;
       float **dist;
       
       dist=declare_float ( T->len, T->len);
    
       for (a=0; a< T->len-1; a++)
	    {
		for ( b=a+1; b< T->len; b++)
		    {
		      
			A=T->ca[a];
			B=T->ca[b];
			dist[a][b]=dist[b][a]=get_atomic_distance ( A,      B);
			
		    }
	    }
       return dist;
   }
float    get_atomic_distance ( Atom *A, Atom*B)
   {
       float dx, dy, dz, d;
       
       if ( !A || !B)
	 {
	   return UNDEFINED;
	 }

       
       dx=A->x - B->x;
       dy=A->y - B->y;
       dz=A->z - B->z;
       

       d=(float) sqrt ( (double) ( dx*dx +dy*dy +dz*dz));
       return d;
   }

char * map_contacts ( char  *file1, char *file2, float T)
{
  
  Ca_trace *ST1, *ST2;
  int *contact_list;
  int a;
  

  ST1=read_ca_trace (file1, "ATOM");
  ST2=read_ca_trace (file2, "ATOM");

  contact_list=identify_contacts (ST1, ST2, T);
  for ( a=0; a<ST1->len; a++)
    {
      ST1->seq[a]=(contact_list[a]==1)?toupper(ST1->seq[a]):tolower(ST1->seq[a]);
    }

  return ST1->seq;
}

float ** print_contacts ( char  *file1, char *file2, float T)
{
  
  Ca_trace *ST1, *ST2;
  float **dist, d;
  int a, b;
  Atom *A, *B;
  char *list;
  int *list1=NULL, *list2=NULL;
  int *cache1, *cache2;
  
  
   
  if ((list=strstr (file1, "_R_"))!=NULL)
      {
	list[0]='\0';
	list+=3;
	list1=string2num_list2(list, "_");
      }
  
  if ((list=strstr (file2, "_R_"))!=NULL)
      {
	list[0]='\0';
	list+=3;
	list2=string2num_list2(list, "_");
      }	

  fprintf ( stdout, "\n#%s (%s) struc2contacts_01", PROGRAM, VERSION);
  fprintf ( stdout, "\nStructure %s vs %s", file1, file2);
  ST1=read_ca_trace (file1, "SEQRES");
  ST2=read_ca_trace (file2, "SEQRES");
  
  
  cache1=vcalloc (ST1->len+1, sizeof (int));
  cache2=vcalloc (ST2->len+1, sizeof (int));
  
  if (list1)for ( a=1; a<list1[0]; a++)cache1[list1[a]]=1;
  else for ( a=1; a<ST1->len; a++)cache1[a]=1;
  
  if (list2)for ( a=1; a<list2[0]; a++)cache2[list2[a]]=1;
  else for ( a=1; a<ST2->len; a++)cache2[a]=1;
  
  
  dist=declare_float (ST1->len+1,ST2->len+1);
  vfree (list1); vfree(list2);
  
  for ( a=0; a< ST1->n_atom; a++)
    {
      A=ST1->structure[a];
      if ( !cache1[A->res_num])continue;
      for ( b=0; b<ST2->n_atom; b++)
	    {
	      
	      B=ST2->structure[b];
	      if( !cache2[B->res_num])continue;
	    
	      d=get_atomic_distance (A,B);
	      
	      if (dist[A->res_num][B->res_num]==0 || dist[A->res_num][B->res_num]>d)dist[A->res_num][B->res_num]=d;
	    }
    }

  for ( a=1; a<=ST1->len; a++)
    {
      A=ST1->ca[a-1];
      if ( !A || !cache1[A->res_num])continue;
      for ( b=1; b<= ST2->len; b++)
	{
	  B=ST2->ca[b-1];
	  if( !B || !cache2[B->res_num])continue;
	  if(dist[a][b]!=0)fprintf ( stdout, "\nResidue %3d [%s] vs %3d [%s] %9.4f Angstrom",A->res_num,A->res,B->res_num,B->res,dist[a][b]);
	}
    }
  fprintf ( stdout, "\n");
  vfree (cache1);vfree(cache2);
  free_float (dist, -1);
  return NULL;
}
  

int * identify_contacts (Ca_trace *ST1,Ca_trace *ST2, float T) 
{
  int a, b;
  float d;
  int *result;
  
  
    
  result=vcalloc ( ST1->len+1, sizeof (int));

  
  for ( a=0; a< ST1->n_atom; a++)
    for ( b=0; b<ST2->n_atom; b++)
  
      {
	
	d=get_atomic_distance (ST1->structure[a], ST2->structure[b]);
	if (d<T)
	  {
	    /*Atom *A, *B;
	    A=ST1->structure[a]; B=ST2->structure[b];
	    fprintf ( stderr, "\n%d %s %s Vs %d %s %s: %f", A->res_num, A->res, A->type, B->res_num, B->res, B->type, d); */
	    result[(ST1->structure[a])->res_num-1]=1;
	  }
      }
  return result;
}

Sequence *seq2contacts ( Sequence *S, float T)
{
  int a;
  Sequence *NS;
  
 
  
  NS=duplicate_sequence (S);
  for ( a=0; a< S->nseq; a++)
    {
      NS->seq[a]=string2contacts ( S->seq[a], S->name[a], S->seq_comment[a], T);
    }
  
  return NS;
}

char *string2contacts (char *seq,char *name, char *comment, float T)
{
  char **nlist;
  char *r;
  char *result;
  int a, b, n;
  char *struc_name;
  static char *struc_file;
  static char *ligand_file;
  Alignment *A;
  char r0, r1;
  
  int l, ln;
  char command[1000];
  /*>seq__struc Ligand1 Chain1 Ligand 2 cahin2
    Chain: index or ANY if unknown
    Ligand: name of pdb file
  */
  
  if ( !struc_file)
    {
      struc_file=vtmpnam  (NULL);
      ligand_file=vtmpnam (NULL);
    }
  

  
  result=vcalloc ( strlen (seq)+1, sizeof (char));
  for ( a=0; a< strlen (seq); a++)result[a]='0';
  
  nlist=string2list (comment);
  if ( !nlist)return result;
  else
    n=atoi(nlist[0]);

  struc_name=strstr(name, "_S_");
  if (!struc_name && !is_pdb_struc (name))
    {
      
      return result;
    }
  else if ( !struc_name && is_pdb_struc (name))
    {
      struc_name=name;
    }
  else
    {
      struc_name+=3;
      if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
	{
	 sprintf (command, "cp %s %s", name, struc_file);
	}
      else
	{
	  sprintf ( command, "extract_from_pdb -infile %s -atom ALL -mode simple -force >%s",name, struc_file);
	}
      my_system (command);
    }
  
  
  
  for ( a=1, ln=1;a<n; ln++)
	{
	  fprintf ( stderr, "\n\tProcess: Structure %s Against Ligand %s T=%.2f Angstrom",struc_name, nlist[a], T);
	 
	  
	  if (check_file_exists ( nlist[a]) && is_simple_pdb_file(nlist[a]))
	      {		
		sprintf (command, "cp %s %s", nlist[a], ligand_file);
		a++;
	      }
	  else if ( check_file_exists ( nlist[a]))
	    {
	      sprintf (command, "extract_from_pdb -infile %s -ligand ALL -atom ALL -mode simple -force >%s", nlist[a], ligand_file);
	      a++;
	    }
	  else
	    {
	      sprintf ( command, "extract_from_pdb -infile %s -chain %s -ligand %s -ligand_only -atom ALL -mode simple -force >%s",struc_name, nlist[a+1],nlist[a], ligand_file);
	      a+=2;	      
	    }
	  my_system (command);
	  
	  if ( T>0)
	    {
	      r=map_contacts (struc_file,ligand_file,T);
	      
	      toggle_case_in_align_two_sequences (KEEP_CASE);
	      A=align_two_sequences (seq,r,"pam250mt", -10, -1, "myers_miller_pair_wise");     
	      toggle_case_in_align_two_sequences (CHANGE_CASE);
	      
	      
	      for ( l=0,b=0; b< A->len_aln; b++)
		{
		  r0=A->seq_al[0][b];r1=A->seq_al[1][b];
		  if (!is_gap(r0))
		    {
		      if (isupper(r1))result[l]=(result[l]!='0')?'9':'0'+ln;
		      l++;
		    }
		}
	      
	      free_aln (A);
	      fprintf ( stderr, " [DONE]");
	    }
	  else if ( T==0)
	    {
	      print_contacts( struc_file,ligand_file,T);
	    }
	  
	}
      fprintf ( stderr, "\n");
      
      return result;
}


char **struclist2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
{
  char *list, **pdb;
  int a;
  char **R, *tmpf;
  
  tmpf=vtmpnam (NULL);
  
  list=strstr ( comment, "_P_")+3;
  if ( !strstr (comment, "_P_"))return NULL;
  else
    {
      pdb=string2list ( strstr ( comment, "_P_")+3);
    }
  
  for (a=1; a<atoi(pdb[0]); a++)
    {
      char com[1000];
      sprintf ( com, "_P_ %s", pdb[a]);
      R=struc2nb (name, seq, com, Threshold, atom, tmpf);
    }
  if (output==NULL|| strstr (output, "fasta"))
    {
      for (a=0; a< strlen (seq); a++)
	if (isupper (R[a][a]))fprintf (stdout, ">%s R=%d T=%.2f %s\n%s\n", name, a+1, Threshold, comment, R[a]);
    }
  else
    {
      FILE *fp;
      char c;
      
      fp=vfopen (tmpf, "r");
      while ( (c=fgetc(fp))!=EOF)fprintf (stdout, "%c", c);
      vfclose (fp);
    }
  return NULL;
}

char **struc2nb (char *name,char *seq, char *comment, float Threshold, char *atom, char *output)
{
  char *struc_file;
  char *struc_name;
  Ca_trace *T;
  Atom *A1, *A2;
  int a, b;
  short **hasch;
  FILE *fp;
  float d;
  char command[10000];
  static char **R;
  
  struc_file=vtmpnam (NULL);
  declare_name (struc_name);
  
  sscanf ( strstr(comment, "_P_"), "_P_ %s", struc_name);
  //struc_name=strstr(name, "_S_");
  
  if (!R)
    {
      int l;
      l=strlen (seq);
      R=declare_char (l+1, l+1);
      for ( a=0; a<l; a++)
	{
	  sprintf ( R[a], "%s", seq);
	  lower_string (R[a]);
	}
    }
  if ( atom) atom=translate_string (atom, "_", " ");
  if ( !struc_name) return NULL;
  
  else 
    {
      //struc_name+=3;
      if ( check_file_exists (struc_name) && is_simple_pdb_file(struc_name))
	{
	  sprintf (command, "cp %s %s", name, struc_file);
	}
      else
	{
	  sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",struc_name,(atom==NULL)?"ALL":atom, struc_file);
	}
  
      my_system (command);
    }
  T=read_ca_trace (struc_file, "ATOM");
  hasch=declare_short (T->len, T->len);
  
  if (!R)
    {
      int l;
      l=strlen (seq);
      R=declare_char (l+1, l+1);
      for ( a=0; a<l; a++)
	{
	  sprintf ( R[a], "%s", seq);
	  lower_string (R[a]);
	}
    }
  for ( a=0; a< T->n_atom; a++)
    for ( b=0; b< T->n_atom; b++)
      {
	A1=T->structure[a];A2=T->structure[b];
	d=get_atomic_distance (A1, A2);

	if ( d<Threshold)
	  {
	    hasch[A1->res_num-1][A2->res_num-1]=1;
	  }
      }
  fp=vfopen (output, "a");
  fprintf ( fp, "#Distance_map_format_01\n#Sequence %s with T= %.2f", struc_name, Threshold);
  for ( a=0; a<T->len; a++)
    {
      int c;
      c=change_residue_coordinate ( T->seq,seq,a);
      
      if ( c!=-1 && seq)
	{
	  char r1, r2;
	  r1=(T->seq)[a];
	  r2=seq[c];
	  r1=tolower(r1);r2=tolower(r2);
	  if ( r1!=r2) continue;
	  R[c][c]=toupper (R[c][c]);
	  fprintf (fp, "\n%s Residue %d ",struc_name,c+1);
	  for ( b=0; b<T->len; b++)
	    {
	      int d;
	      char r3, r4;
	      r3=(T->seq)[a];r4=seq[c];
	      r3=tolower(r3);r4=tolower(r4);
	      if ( r3!=r4) continue;
	      
	      d=change_residue_coordinate (T->seq,seq,b);
	      
	      if ( hasch[a][b] && d!=-1)
		{
		  fprintf (fp, "%d ",d+1);
		  R[c][d]=toupper (R[c][d]);
		}
	    }
	  fprintf (fp, ";");
	}
    }
  free_short (hasch, -1);
  fprintf (fp, "\n");
  return R;
}

short **seq2nb (char *seq, char *pdb, float Threshold, char *atom)
{
  char *struc_file;
  char *struc_name;
  Ca_trace *T;
  Atom *A1, *A2;
  int a, b;
  short **hasch;
  short **result; //Contains the result for every residue of seq
  float d;
  char command[10000];

  // get a clean pdb file
  struc_file=vtmpnam (NULL);
  if ( check_file_exists (struc_file) && is_simple_pdb_file(struc_name))
    {
      sprintf (command, "cp %s %s", pdb, struc_file);
    }
  else
    {
      sprintf ( command, "extract_from_pdb -infile %s -atom %s -mode simple -force >%s",pdb,(atom==NULL)?"ALL":atom, struc_file);
    }
  my_system (command);
  
  //Read and hasch the PDB file
  T=read_ca_trace (struc_file, "ATOM");
  hasch=declare_short (T->len, T->len);
  result=declare_short (strlen (seq)+1, strlen (seq)+1);
  for ( a=0; a< T->n_atom; a++)
    for ( b=0; b< T->n_atom; b++)
      {
	A1=T->structure[a];A2=T->structure[b];
	d=get_atomic_distance (A1, A2);

	if ( d<Threshold)
	  {
	    hasch[A1->res_num-1][A2->res_num-1]=1;
	  }
      }
  
  for ( a=0; a<T->len; a++)
    {
      int c;
      c=change_residue_coordinate ( T->seq,seq,a);
      if ( c!=-1)
	{
	  for ( b=0; b<T->len; b++)
	    {
	      int d;
	      d=change_residue_coordinate (T->seq,seq,b);
	      
	      if ( hasch[a][b] && d!=-1)
		{
		  result[c][result[c][0]++]=d;
		}
	    }
	}
    }
  free_short (hasch, -1);
  return result;
}
/*********************************COPYRIGHT NOTICE**********************************/
/*© Centro de Regulacio Genomica */
/*and */
/*Cedric Notredame */
/*Tue Oct 27 10:12:26 WEST 2009. */
/*All rights reserved.*/
/*This file is part of T-COFFEE.*/
/**/
/*    T-COFFEE 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.*/
/**/
/*    T-COFFEE is distributed in the hope that it will be useful,*/
/*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
/*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
/*    GNU General Public License for more details.*/
/**/
/*    You should have received a copy of the GNU General Public License*/
/*    along with Foobar; if not, write to the Free Software*/
/*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
/*...............................................                                                                                      |*/
/*  If you need some more information*/
/*  cedric.notredame@europe.com*/
/*...............................................                                                                                                                                     |*/
/**/
/**/
/*	*/
/*********************************COPYRIGHT NOTICE**********************************/