/*
  Program RMS   written by Scott M. Le Grand
  Even more modified and completely changed in 1999 by Arne Elofsson 
  Modified by Arne Elofsson in 1994
  Copyright 1992, Scott M. Le Grand
  */

#include <stdio.h>
#include <math.h> 
#ifdef  NOZLIB
#else
#include <zlib.h> 
#endif
#include <sys/types.h>
#include <sys/times.h>
#include <sys/param.h>
#include <sys/time.h> 

#define	MAXATMS		5000		/* Maximum allowable atoms */
#define PI		3.14159265	/* Useful constant */

#define TRUE		1		/* Boolean definitions */
#define FALSE		0

typedef struct {
    struct {
	double x,y,z;		/* Atomic coordinates */
	double rms;		/* RMS deviation */
	char residue[8];	/* PDB info for output */
	char name[8];
	int number;
	int resnum;
	int selected;
    } atm[MAXATMS];
    double xcen,ycen,zcen;
    int	atoms;			/* Current # of atoms */
    char	filename[80];		/* filename to read molecule from */
} molecule;

molecule	m[2];		/* Molecules to be compared */

double		LG_pvalue();
double		Levitt_Gerstein();
double		superimpose_molecules();
int		b_flag;		
int		align_flag;		
int		select1_flag;		
int		select2_flag;		
int		verbose_flag;		/* Verbose output flag */
int             output_flag;            /* Output flag */
int             fraction_flag;            /* Output flag */
char            output_name[80];            /* Output name */
int             output_flag_ca;            /* Output flag */
char            output_name_ca[80];            /* Output name */
int             read_molecules_ca();
double		rms;			/* final RMS error */
/*int             temp;*/

main(argc,argv)		/* Main routine */
    int argc;
    char *argv[];
{
 molecule      m0,m1;
 int	files;	/* Number of files parsed */
 int	i,j,k;	/* Counter variable */
 int	L=4;	// Minimum length of fragment for FindMaxSub matching
 int    minatoms=25; /* Variable for min no of atoms for P-value*/
 double	s[3][3];	/* Final transformation matrix */
 double          maxrmsd,minrmsd;
 double          fraction=1;
 double          numfrac=0;
 double          maxpvalue=1,maxpvalue1=1.;
 double          score,maxrms,maxrms1,maxscore,maxscore1;
 double          pvalue;
 int             maxatoms,maxpos,minpos,atoms;
 int             maxsub_15=0;
 int             maxsub_20=0;
 int             maxsub_25=0;
 int             maxsub_30=0;
 int             maxsub_35=0;
 int             maxselect1_0[MAXATMS],maxselect1_1[MAXATMS];
 int             maxselect2_0[MAXATMS],maxselect2_1[MAXATMS];
 double          bestsizescore[MAXATMS];
 int             maxatoms1;
 /* Initialize */
 files=0;
 verbose_flag=FALSE;
 output_flag=FALSE;
 select2_flag=FALSE;
 select1_flag=FALSE;
      
 /* Parse command line for PDB files and options */
 i=1;
/*    temp=1; */
 while (i<argc)
  {
    if (strcmp(argv[i],"-v")==0)
      {
       verbose_flag=TRUE;
       select2_flag=TRUE;
       select1_flag=TRUE;
      }
    else if (strcmp(argv[i],"-f")==0)
      {
	numfrac=0.2;
	fraction_flag=TRUE;
      }
     else if (strcmp(argv[i],"-1")==0)
      {
	select1_flag=TRUE;
      }
    else if (strcmp(argv[i],"-2")==0)
      {
	select2_flag=TRUE;
      }
    else if (strcmp(argv[i],"-o")==0)
      {
	align_flag=TRUE;
      }
    else if (strcmp(argv[i],"-b")==0)
      {
	b_flag=TRUE;
      }
    else if (strcmp(argv[i],"-a")==0)
      {
	numfrac=0.;
	fraction_flag=TRUE;
	select2_flag=TRUE;
	select1_flag=TRUE;
      }
    else
      {
       strcpy(m[files].filename,argv[i]);
       files++;
     }
   i++;
  }
 if (files==2)	/* Valid input line? */
   {
     if (read_molecules_ca()==0)
       {
	 check_molecules(&m[0],&m[1]); // this deletes some files
	 //printf("testing %d %d\n",m[0].atoms,m[1].atoms);
	 atoms=m[0].atoms;


	 /*
	 center_molecule(&m[0]);
	 center_molecule(&m[1]);
	 rms=superimpose_molecules(&m[0],&m[1],s);
	 printf("test %f\n",rms);
	 exit();
	 */
	 /* ************************************************************************
	    First we'll try to do this the otherway, i.e. adding atoms
	  ************************************************************************ */
	 for (i=0;i<m[0].atoms;i++){bestsizescore[i]=0;}
	 for (i=0;i<m[0].atoms-L;i++){
	   for (j=0;j<m[0].atoms;j++){ // deselect all atoms
	     m[0].atm[j].selected=FALSE;
	     m[1].atm[j].selected=FALSE;
	   }
	   for (j=i;j<i+L;j++){ // select some
	     m[0].atm[j].selected=TRUE;
	     m[1].atm[j].selected=TRUE;
	   }
	   j=L;
	   center_molecule(&m[0]);
	   center_molecule(&m[1]);
	   rms=superimpose_molecules(&m[0],&m[1],s);	   
	   while (j<m[0].atoms){
	     minrmsd=9999.0;
	     for (k=0;k<m[0].atoms;k++) {
	       if (m[0].atm[k].rms < minrmsd && ! m[0].atm[k].selected ){
		 minpos=k;
		 minrmsd=m[0].atm[k].rms;
	       }
	     }
	     j++;
	     //	     printf("Selecting %d %f\n",minpos,minrmsd);
	     m[0].atm[minpos].selected=TRUE;
	     m[1].atm[minpos].selected=TRUE;
	     //	     printf("testing2a %f\n",m[1].atm[i].x);
	     center_molecule(&m[0]);
	     center_molecule(&m[1]);
	     rms=superimpose_molecules(&m[0],&m[1],s);
	     score=Levitt_Gerstein(&m[0],&m[1]);
	     pvalue=LG_pvalue(j,score);
	     if (b_flag){
	       if (score > bestsizescore[j]){bestsizescore[j]=score;}
	     }

	     if (rms<1.5 && maxsub_15<j){maxsub_15=j;}
	     if (rms<2.0 && maxsub_20<j){maxsub_20=j;}
	     if (rms<2.5 && maxsub_25<j){maxsub_25=j;}
	     if (rms<3.0 && maxsub_30<j){maxsub_30=j;}
	     if (rms<3.5 && maxsub_35<j){maxsub_35=j;}


	     if (fraction_flag){
	       if ((double)j/atoms <= fraction ){
		 printf("FRAC1:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
			(double)j/atoms,j,rms,score,
			pvalue);
		 fraction-=numfrac;
		 //printf("TEST %f %f\n",fraction,numfrac);
	       } 
	     }

	     if (verbose_flag) output_results();
	     if ((pvalue <= maxpvalue) && (j>minatoms)){
	       if (fraction_flag){
		 printf("GOOD1:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
			(double)j/atoms,j,rms,score,
			pvalue);
	       }
	       //printf ("TEST:\t%d\t%f\t%f\t%f\n",m[0].atoms,pvalue,rms,score);
	       maxpvalue=pvalue;
	       maxatoms=j;
	       maxrms=rms;
	       maxscore=score;
	       for (k=0;k<m[0].atoms;k++){maxselect1_0[k]=m[0].atm[k].selected;}
	       for (k=0;k<m[1].atoms;k++){maxselect1_1[k]=m[1].atm[k].selected;}
	     }
	     //	     printf("testing5 %f\n",m[1].atm[i].x);
	   }
	 }
	 //	 printf("BEST1:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
	 // (double)maxatoms/(double)atoms,maxatoms,maxrms,maxscore,maxpvalue);
	 
	 maxatoms1=maxatoms;
	 maxrms1=maxrms;
	 maxscore1=maxscore;
	 maxpvalue1=maxpvalue;


	 /* ************************************************************************
	  Then we try to delete atoms
	  ************************************************************************ */
	 maxpvalue=1.1;
	 fraction=1;
	 for (i=0;i<m[0].atoms;i++){ // select all atoms
	   m[0].atm[i].selected=TRUE;
	   m[1].atm[i].selected=TRUE;
	 }
	 atoms=m[0].atoms;

	 while (atoms>5){
	   center_molecule(&m[0]);
	   center_molecule(&m[1]);
	   rms=superimpose_molecules(&m[0],&m[1],s);
	   score=Levitt_Gerstein(&m[0],&m[1]);
	   pvalue=LG_pvalue(atoms,score);
	   if (b_flag){
	     if (score > bestsizescore[atoms]){bestsizescore[atoms]=score;}
	   }

	   if (rms<1.5 && maxsub_15<atoms){maxsub_15=atoms;}
	   if (rms<2.0 && maxsub_20<atoms){maxsub_20=atoms;}
	   if (rms<2.5 && maxsub_25<atoms){maxsub_25=atoms;}
	   if (rms<3.0 && maxsub_30<atoms){maxsub_30=atoms;}
	   if (rms<3.5 && maxsub_35<atoms){maxsub_35=atoms;}

	   if (m[0].atoms == atoms  ){
	     if (fraction_flag){
	       printf("ALL:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
		      (double)atoms/m[0].atoms,atoms,rms,score, pvalue);
	     }
	     maxpvalue=pvalue;
	     maxatoms=atoms;
	     maxrms=rms;
	     maxscore=score;	    
	     for (k=0;k<m[0].atoms;k++){maxselect2_0[k]=m[0].atm[k].selected;}
	     for (k=0;k<m[1].atoms;k++){maxselect2_1[k]=m[1].atm[k].selected;}
	   }

	   if (fraction_flag){
	     if ((double)atoms/m[0].atoms <= fraction ){
	       printf("FRAC2:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",	
		      (double)atoms/m[0].atoms,atoms,rms,score, pvalue);
	       fraction-=numfrac;
	       //	     printf("TEST %f %f\n",fraction,numfrac);
	     } 
	   }

	   if (verbose_flag) output_results();
	   if ((pvalue < maxpvalue) && (atoms>minatoms)){
	     //	     printf ("TEST:\t%d\t%f\t%f\t%f\n",m[0].atoms,pvalue,rms,score);
	     if (fraction_flag){
	       printf("GOOD2:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",	
		      (double)atoms/m[0].atoms,atoms,rms,score, pvalue);
	     }
	     fraction-=numfrac;
	     maxpvalue=pvalue;
	     maxatoms=atoms;
	     maxrms=rms;
	     maxscore=score;
	     for (k=0;k<m[0].atoms;k++){maxselect2_0[k]=m[0].atm[k].selected;}
	     for (k=0;k<m[1].atoms;k++){maxselect2_1[k]=m[1].atm[k].selected;}
	   }
	   maxrmsd=-0.1;
	   maxpos=-1;
	   for (i=0;i<m[0].atoms;i++)
	     {
	       //printf("TEST-y:\t %d %d %f %f\n",i,m[0].atm[i].selected,m[0].atm[i].rms,maxrmsd);
	       if (m[0].atm[i].selected && m[0].atm[i].rms > maxrmsd){
		 maxpos=i;
		 maxrmsd=m[0].atm[i].rms;
	       }
	     }
	   //	   delete_atom(&m[0],maxpos);
	   // delete_atom(&m[1],maxpos);
	   //printf("Unselecting:\t%d\t%f\n",maxpos,maxrmsd);
	   if(maxpos<0){printf ("ERROR2, not selecting any atom %d\n",i);}
	   else{
	     m[0].atm[maxpos].selected=FALSE;
	     m[1].atm[maxpos].selected=FALSE;
	   }
	   atoms--;
	 }
	 //printf("BEST2:\t%.2f\t%d\t%.1f\t%.1f\t%e\t%e\n",
	 //(double)maxatoms/(double)m[0].atoms,maxatoms,maxrms,maxscore,maxpvalue,maxpvalue1);
	 if (maxpvalue >= maxpvalue1 || select1_flag){
	   if (maxatoms1<minatoms){maxatoms1=maxrms1=maxscore1=0;maxpvalue1=1;}
	   printf("SCORE-1:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
	   (double)maxatoms1/(double)m[0].atoms,maxatoms1,maxrms1,maxscore1,maxpvalue1);
	   if (align_flag){
	     for (i=0;i<m[0].atoms;i++){m[0].atm[i].selected=maxselect1_0[i];}
	     for (i=0;i<m[1].atoms;i++){m[1].atm[i].selected=maxselect1_1[i];}
	     center_molecule(&m[0]);
	     center_molecule(&m[1]);
	     rms=superimpose_molecules(&m[0],&m[1],s);	   
	     output_results();
	   }
	 }
	 if (maxpvalue < maxpvalue1 || select2_flag){
	   if (maxatoms<minatoms){maxatoms=maxrms=maxscore=0;maxpvalue=1;}
	   printf("SCORE-2:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
		(double)maxatoms/(double)m[0].atoms,maxatoms,maxrms,maxscore,maxpvalue);
	   if (align_flag){
	     for (i=0;i<m[0].atoms;i++){m[0].atm[i].selected=maxselect2_0[i];}
	     for (i=0;i<m[1].atoms;i++){m[1].atm[i].selected=maxselect2_1[i];}
	     center_molecule(&m[0]);
	     center_molecule(&m[1]);
	     rms=superimpose_molecules(&m[0],&m[1],s);	   
	     output_results();
	   }
	 }
	 //MAXSUB:   \t1.5\t2.0\t2.5\t3.0\t3.5\n
	 printf("MAXSUB:   \t%d\t%d\t%d\t%d\t%d\n",
		maxsub_15,maxsub_20,maxsub_25,maxsub_25,maxsub_30);
		
	 
	 if (b_flag){
	   for (i=0;i<m[0].atoms;i++){
	     printf("SCORE-TEST\t%d\t%f\t%e\n",i,bestsizescore[i],LG_pvalue(i,bestsizescore[i]));
	   }
	 }
       }
     else
       printf("Usage: rms [ -a -f -v -1 -2 ] file1 file2\n");
   }
}



int read_molecules_ca(temp)	/* Reads in molecules to be superimposed */
    int temp;
{
  int	i,j,k,atoms;	/* Counter variables */
  char	buff[512];	/* Input string */
  char	junk[512];	/* Throwaway string */
  char	line_flag[512];	/* PDB file line mode */
  char	residue[8];	/* PDB atom info for output */
  char	name[8];
  int	resnum;
  int	number;
  double	x,y,z;		/* Temporary coordinates values */
  FILE *fp;
    
  for (i=0;i<2;i++)
    {
#ifdef  NOZLIB
      fp=fopen(m[i].filename,"r");	/* Does file exist? */
#else
      fp=gzopen(m[i].filename,"r");	/* Does file exist? */
#endif
      if (fp!=NULL)	/* If yes, read in coordinates */
	{
	  /* Initialize things */
	  m[i].xcen=m[i].ycen=m[i].zcen=0;
	  atoms=0;
#ifdef NOZLIB
	  while(fgets(buff,255,fp)!=NULL)
#else
	  while(gzgets(fp,buff,255)!=NULL)
#endif
	    {
	      sscanf(buff,"%s %d %s %s",line_flag,&number,name,residue);
	      if (strcmp("ATOM",line_flag)==0)	/* Is it an ATOM entry? */
		{
		  if (strcmp("CA",name)==0)	   /*  Is it an CA atom ? */
		    {
		      sscanf(&buff[22],"%d %lf %lf %lf",&resnum,&x,&y,&z);
		      //		      printf("test %d\n",resnum);
		      m[i].atm[atoms].x=x;
		      m[i].atm[atoms].y=y;
		      m[i].atm[atoms].z=z;
		      m[i].atm[atoms].resnum=resnum;
		      m[i].atm[atoms].number=number;
		      m[i].atm[atoms].selected=TRUE;
		      strcpy(m[i].atm[atoms].name,name);
		      strcpy(m[i].atm[atoms].residue,residue);
		      m[i].xcen+=x;
		      m[i].ycen+=y;
		      m[i].zcen+=z;
		      atoms++;
		    }
		}
	    }
	  m[i].atoms=atoms;
#ifdef NOZLIB
	  close(fp);
#else
	  gzclose(fp);
#endif
	  //	  if (atoms!=m[0].atoms)		/* Are file sizes indentical? */
	  //{
	  //  printf("Inconsistent number of atoms in file %s\n",m[i].filename);
	  //  return(1);
	  //}
	
	}
      else
	{
	  printf("Couldn't open file %s\n",m[i].filename);
	  exit();
	}
    }
  return(0);
}
int center_molecule(m)	/* Reads in molecules to be superimposed */
    molecule *m;
{
  int	i;	/* Counter variables */
  int   natoms;  // Number of selected atoms
  double	xcen,ycen,zcen;		/* Temporary coordinates values */

  natoms=0;
  xcen=ycen=zcen=0;
  for (i=0;i<m->atoms;i++){
    if (m->atm[i].selected){
      xcen+=m->atm[i].x;
      ycen+=m->atm[i].y;
      zcen+=m->atm[i].z;
      natoms++;
    }
  }
  /* Now center molecule */
  xcen/=(double)natoms;
  ycen/=(double)natoms;
  zcen/=(double)natoms;
  //printf("TEST natoms: %d %f %f %f \n",natoms,xcen,ycen,zcen);
  for (i=0;i<m->atoms;i++)
    {
      m->atm[i].x-=xcen;
      m->atm[i].y-=ycen;
      m->atm[i].z-=zcen;
    }
  return(0);
}
multiply_matrix(a,b,c)  /* computes C=AB */
    double a[3][3],b[3][3],c[3][3];
{
    int i,j,k;
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    {
		c[i][j]=0.0;
		for (k=0;k<3;k++)
		    c[i][j]+=a[i][k]*b[k][j];
	    }
}
copy_matrix(f,t) /* copy matrix f into matrix t */
    double f[3][3],t[3][3];
{
    int i,j;
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    t[i][j]=f[i][j];
}
transpose_matrix(m) /* Transpose a 3x3 matrix */
    double m[3][3];
{
    double dummy;
    dummy=m[0][1]; m[0][1]=m[1][0]; m[1][0]=dummy;
    dummy=m[0][2]; m[0][2]=m[2][0]; m[2][0]=dummy;
    dummy=m[1][2]; m[1][2]=m[2][1]; m[2][1]=dummy;
}
delete_atom(m1,num)	
    molecule 	*m1;		
    int         num;
{
  int          i,j,k;
  for (i=num;i<m1->atoms;i++)
    {
      j=i+1;
      strcpy(m1->atm[i].name,m1->atm[j].name);
      strcpy(m1->atm[i].residue,m1->atm[j].residue);
      m1->atm[i].x=m1->atm[j].x;
      m1->atm[i].y=m1->atm[j].y;
      m1->atm[i].z=m1->atm[j].z;
      /*      m1->atm[i].rms=m1->atm[j].rms; */
      m1->atm[i].number=m1->atm[j].number;
      m1->atm[i].resnum=m1->atm[j].resnum;

    }
  m1->atoms--;
}

check_molecules(m1,m2)	
    molecule 	*m1,*m2;		
{
  /* this will delete all atoms that are not in both molecules */
  /* Now we should extract the part of the molecules that exists in both */
  int          i,j,minlength;
  int          found;

  minlength=m1->atoms;
  if (minlength>m1->atoms) minlength=m2->atoms;
  i=0;
  while (i<=m1->atoms &&  i<=m2->atoms ){
    while ( m1->atm[i].resnum < m2->atm[i].resnum && i< m1->atoms ){
      //      printf("Deleting m1: %d %d %d\n",i,m1->atm[i].resnum,m2->atm[i].resnum);
      delete_atom(m1,i);	
    }
    while (m1->atm[i].resnum > m2->atm[i].resnum &&  i< m2->atoms ){
      //	printf("Deleting m2: %d %d %d\n",i,m1->atm[i].resnum,m2->atm[i].resnum);
      delete_atom(m2,i);	
    }
    i++;
  }
  while (m1->atoms < m2->atoms){    delete_atom(m2,m2->atoms)  ;}
  while (m1->atoms > m2->atoms){    delete_atom(m1,m1->atoms)  ;}
}


double Levitt_Gerstein(m1,m2)	
    molecule 	*m1,*m2;		
{
  int          i,j,last1,last2;
  int          numgap=0;
  double       sum=0.;
  //  double       d0=25.; // 5**2
  double       d0=5.;  // 2.24**2
  double       M=20.;
  last1=m1->atm[0].resnum-1;
  last2=m2->atm[0].resnum-1;
  j=0;
  for (i=0;i<m1->atoms;i++)
    {
      if (m1->atm[i].selected){
	last1++;
	last2++;
	if (m1->atm[i].resnum > last1 || m2->atm[i].resnum > last2 ){
	  numgap++;
	}
	sum+=(1/(1+m1->atm[i].rms/d0));
	//	printf("Test>\t%d\t%f\t%f\n",numgap,sum,m1->atm[i].rms);
	last1=m1->atm[i].resnum;
	last2=m2->atm[i].resnum;
	j++;
	//	printf("Test> %d  \t%d\t%d\t%f\n",i,m1->atm[i].selected,sum);
      }
    }
  //printf("TEST-LG>\t%d\t%d\t%f\t%f\n",j,numgap,sum,M*(sum-numgap/2));
  return(M*(sum-numgap/2));
}

double LG_pvalue(N,MS)
     int N;
     double MS;
{
  double min_SD=1.e-20;
  double ln,Mean_MS,SD_MS,Z,expZ;
  double ln60=log(120.);
  double b=2.0*ln60*18.411424-4.501719;
  double a=(ln60*ln60)*18.411424+ln60*(-4.501719)+2.637163-b*ln60;
  if (N<6) {
    expZ=1.;
    return(expZ);
  }
  ln=log((double)N);
  if (N<120){
    Mean_MS=ln*ln*18.411424+ln*(-4.501719)+ 2.637163;
    SD_MS=ln*21.351857 -37.521269;
  }else{
    Mean_MS=ln*b+a;
    SD_MS=ln60*21.351857 -37.521269;
  }
  //  if (SD_MS<min_SD) {SD_MS=min_SD;}
  Z=(MS-Mean_MS)/SD_MS;

  expZ=exp((double)-1.*Z);
  //printf("TEST-P> %d\t%f\t%f\t%e\t%e\t%e\t%e\n",N,MS,Mean_MS,SD_MS,Z,expZ,1-exp(-1*expZ));
  if (Z>20){
    return(expZ);
    //  }else if (Z<-100){
    //expZ=1.;
    //return(expZ);
  }else{
    return(1-exp(-expZ));
  }
}

double superimpose_molecules(m1,m2,s)	/* Find RMS superimposition of m1 on m2 */
    molecule 	*m1,*m2;		/* Molecules to be superimposed */
    double 		s[3][3];	/* Final transformation matrix */
{
    int 		i,j,k;		/* Counter variables */
    int                 natoms=0;
    double		u[3][3];	/* direct product matrix */
    double 		t[3][3];	/* Temporary storage matrix */
    double		ma[3][3];	/* x axis rotation matrix */
    double		mb[3][3];	/* y axis rotation matrix */
    double		mg[3][3];	/* z axis rotation matrix */
    double 		*d1,*d2;	/* usefule pointers */
    double 		error;		/* Final superimposition error */
    double 		error2;		/* Final superimposition error */
    double		alpha=0.0; 	/* Angle of rotation around x axis */
    double		beta=0.0;	/* Angle of rotation around y axis */
    double		gamma=0.0;	/* Angle of rotation around z axis */
    double		dist,dist2,x,y,z;		/* Temporary coordinate variables */
    for (i=0;i<3;i++)  /* Initialize matrices */
	for (j=0;j<3;j++)
	    s[i][j]=u[i][j]=0.0;
    s[0][0]=s[1][1]=s[2][2]=1.0;  /* Initialize S matrix to I */
    for (i=0;i<3;i++)  /* Initialize rotation matrices to I */
	for (j=0;j<3;j++)
	    ma[i][j]=mb[i][j]=mg[i][j]=s[i][j];
    for (i=0;i<m1->atoms;i++)  /* Construct U matrix */
	{
	  if (m1->atm[i].selected){
	    d1= &(m1->atm[i].x);
	    for (j=0;j<3;j++)
	      {
		d2= &(m2->atm[i].x);
		for (k=0;k<3;k++)
		  {
		    u[j][k]+=(*d1)*(*d2);
		    d2++;
		  }
		d1++;
	      }
	  }
	}
	  
    do
      {
	error=0.0;
	/* Calculate x axis rotation */
	alpha=atan((u[2][1]-u[1][2])/(u[1][1]+u[2][2]));
	/* Insure we are heading for a minimum, not a maximum */
	if (cos(alpha)*(u[1][1]+u[2][2])+sin(alpha)*(u[2][1]-u[1][2])<0.0)
	  alpha+=PI;
	ma[1][1]=ma[2][2]=cos(alpha);
	ma[2][1]=sin(alpha); ma[1][2]= -ma[2][1];
	transpose_matrix(ma);
	multiply_matrix(u,ma,t);
	transpose_matrix(ma);
	copy_matrix(t,u);
	multiply_matrix(ma,s,t);
	copy_matrix(t,s);
	/* Calculate y axis rotation */
	beta=atan((u[0][2]-u[2][0])/(u[0][0]+u[2][2]));
	/* Insure we are heading for a minimum, not a maximum */
	if (cos(beta)*(u[0][0]+u[2][2])+sin(beta)*(u[0][2]-u[2][0])<0.0)
	  beta+=PI;
	mb[0][0]=mb[2][2]=cos(beta);
	mb[0][2]=sin(beta); mb[2][0]= -mb[0][2];
	transpose_matrix(mb); 
	multiply_matrix(u,mb,t);
	transpose_matrix(mb);
	copy_matrix(t,u);
	multiply_matrix(mb,s,t);
	copy_matrix(t,s); 
	/* Calculate z axis rotation */
	gamma=atan((u[1][0]-u[0][1])/(u[0][0]+u[1][1]));
	/* Insure we are heading for a minimum, not a maximum */
	if (cos(gamma)*(u[0][0]+u[1][1])+sin(gamma)*(u[1][0]-u[0][1])<0.0)
	  gamma+=PI;
	mg[0][0]=mg[1][1]=cos(gamma);
	mg[1][0]=sin(gamma); mg[0][1]= -mg[1][0];
	transpose_matrix(mg);
	multiply_matrix(u,mg,t);
	transpose_matrix(mg);
	copy_matrix(t,u);
	multiply_matrix(mg,s,t);
	copy_matrix(t,s);
	error=fabs(alpha)+fabs(beta)+fabs(gamma); 
      }
    while (error>0.0001);	/* Is error low enough to stop? */
    /* Now calculate final RMS superimposition */
    error=0.0;
    error2=0.0;
    natoms=0;
    /*    printf("testing7b %f %f\n",m2->atm[1].x,error);
    printf ("s1:\t%f\t%f\t%f\n",s[0][0],s[0][1],s[0][2]);
    printf ("s2:\t%f\t%f\t%f\n",s[1][0],s[1][1],s[1][2]);
    printf ("s3:\t%f\t%f\t%f\n",s[2][0],s[2][1],s[2][2]);*/
    if (! (s[0][0]>0 || s[0][0]<0))  {
      //      printf ("bar\n");
      for (i=0;i<3;i++)  /* Initialize matrices */
	for (j=0;j<3;j++)
	  s[i][j]=u[i][j]=0.0;
      s[0][0]=s[1][1]=s[2][2]=1.0;  /* Initialize S matrix to I */
    }
    for (i=0;i<m2->atoms;i++)
      {

	dist=m2->atm[i].x*m2->atm[i].x+m2->atm[i].y*m2->atm[i].y+m2->atm[i].z*m2->atm[i].z;
	//	printf("TEST1-m2\t %d\t%6.3f %6.3f %6.3f\t%e\n",i,m2->atm[i].x,m2->atm[i].y,m2->atm[i].z,dist);
	x=s[0][0]*m2->atm[i].x+s[0][1]*m2->atm[i].y+s[0][2]*m2->atm[i].z;
	y=s[1][0]*m2->atm[i].x+s[1][1]*m2->atm[i].y+s[1][2]*m2->atm[i].z;
	z=s[2][0]*m2->atm[i].x+s[2][1]*m2->atm[i].y+s[2][2]*m2->atm[i].z;
	m2->atm[i].x=x;
	m2->atm[i].y=y;
	m2->atm[i].z=z;
	//dist2=m2->atm[i].x*m2->atm[i].x+m2->atm[i].y*m2->atm[i].y+m2->atm[i].z*m2->atm[i].z;
	//printf("TEST2-m2\t %d\t%6.3f %6.3f %6.3f\t%e\t%e\n",i,m2->atm[i].x,m2->atm[i].y,m2->atm[i].z,dist2,dist2-dist);
	x=m1->atm[i].x-x;
	y=m1->atm[i].y-y;
	z=m1->atm[i].z-z;
	m2->atm[i].rms=m1->atm[i].rms=x*x+y*y+z*z;
	error+=m1->atm[i].rms;
	if (m1->atm[i].selected){
	  error2+=m1->atm[i].rms;
	  natoms++;
	}
	/*printf ("TEST1: %f %f %f\n",m1->atm[i].x,m1->atm[i].y,m1->atm[i].z);
		  printf ("TEST2: %f %f %f\n",m2->atm[i].x,m2->atm[i].y,m2->atm[i].z);
		  printf ("TEST3: %d %f %f\n",i,error,m1->atm[i].rms); */
      }
    //    printf("testing8 %f\n",m2->atm[1].x);
    error/=(double)(m1->atoms);
    error2/=(double)(natoms);
    return(sqrt(error2));
}
output_results()
{
  int	i,j;		/* Counter variable */
  double	sum;		/* Average RMS counter by residue */
  double	resatoms;	/* Number of atoms in each residue */
  char	name[8];	/* Temporary residue name variable */
  printf("Overall average RMS deviation is %lf Angstroms.\n",rms);
  /* First output RMS deviations by atom */
  printf("\n");
  printf("RMS deviations by atom\n");
  printf("Atom \t Name \t  residue \t RMS deviation\tAtom \t Name \t  residue \t RMS deviation\n");
  for (i=0;i<m[0].atoms;i++){
    /*	  printf ("Test1: %f %f %f\n",m[0].atm[i].x,m[0].atm[i].y,m[0].atm[i].z);
	  printf ("Test2: %f %f %f\n",m[1].atm[i].x,m[1].atm[i].y,m[1].atm[i].z); */
    if ( m[0].atm[i].selected){
      printf("%d \t %s \t %3s %4d \t %8.4lf\t%d \t %s \t %3s %4d \t %8.4lf\n",
	     m[0].atm[i].number,m[0].atm[i].name,
	     m[0].atm[i].residue,m[0].atm[i].resnum,
	     sqrt(m[0].atm[i].rms),
	     m[1].atm[i].number,m[1].atm[i].name,
	     m[1].atm[i].residue,m[1].atm[i].resnum,
	     sqrt(m[1].atm[i].rms));
    }
  }
}
output_file(filename)
    char filename[80];
{
    int i;
    FILE *fp;
    fp=fopen(filename,"w");
    if (fp!=NULL)
	{
	    for (i=0;i<m[1].atoms;i++)
	      {
		fprintf(fp,"ATOM   %4d  %-4s%3s  %4d     %7.3lf %7.3lf %7.3lf  1.0 %5.2lf      RMSD\n",
			m[1].atm[i].number,m[1].atm[i].name,
			m[1].atm[i].residue,m[1].atm[i].resnum,
			m[1].atm[i].x,m[1].atm[i].y,m[1].atm[i].z,
			sqrt(m[0].atm[i].rms));
	      }
	    fprintf(fp,"END \n");
/*	    printf("test %d  %d \n",m[1].atoms,m[0].atoms);*/
	    for (i=0;i<m[0].atoms;i++)
	      {
		fprintf(fp,"ATOM   %4d  %-4s%3s  %4d     %7.3lf %7.3lf %7.3lf  1.0 %5.2lf      RMSD\n",
			m[0].atm[i].number,m[0].atm[i].name,
			m[0].atm[i].residue,m[0].atm[i].resnum,
			m[0].atm[i].x,m[0].atm[i].y,m[0].atm[i].z,
			sqrt(m[0].atm[i].rms));
	      }
	    fprintf(fp,"END \n");
	}
    close(fp);
}