/*
  Program LGSCORE2   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 and Arne Elofsson (1999)
  The program is available under the Gnu public license (see www.fsf.org)

 This program 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; version 2 of the License. With the
 exception that if you use this program in any scientific work you have
 to explicitly state that you have used LGSCORE and cite the relevant
 publication (dependent on what you have used LGSCORE for). My
 publicationlist can be found at http://www.sbc.su.se/~arne/papers/.
 This program 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 this program (in 
 the file gpl.txt); if not, write to the Free Software Foundation, 
 Inc., 59 Temple Place - Suite 330, Boston, MA 02111-13

  */

#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		1000		/* Maximum allowable atoms */
#define PI		3.14159265	/* Useful constant */
#define TRUE		1		/* Boolean definitions */
#define FALSE		0
#define d0              5.00000  // 2.24**2
#define M               20.00000
#define penopen         -10.00000
#define max(A,B) ((A) > (B) ? (A) : (B))
#define min(A,B) ((A) < (B) ? (A) : (B))

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[800];		/* filename to read molecule from */
} molecule;

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

double		LG_pvalue();
double		LG_pvalueF();
double		Levitt_Gerstein();
double		superimpose_molecules();
int	shortest;			/* Current # of atoms */
int		verbose_flag;		/* Verbose output flag */
int		select1_flag;
int		select2_flag;
int		b_flag;
int		align_flag;
int             output_flag;            /* Output flag */
int             out_flag;            /* Output flag */
int             fraction_flag;            /* Output flag */
char            output_name[800];            /* Output name */
int             output_flag_ca;            /* Output flag */
char            output_name_ca[800];            /* Output name */
char            outfile[800];            /* Output name */
int             read_molecules_ca();
int             nalign;    // number of aligned residues
double		rms;			/* final RMS error */
double          ss;
double          dyntrace();
double          sdist();
//double          max();
//double          min();
int 		length,ma[2];
/*int             temp;*/
 double          cutoffrmsd=2.5;
 double          cutoffrms=6;
 double          cutoff;
 int             cut;
 int             count;
 int             maxcount=10;

main(argc,argv)		/* Main routine */
    int argc;
    char *argv[];
{
 int	files;	/* Number of files parsed */
 int	a,b,i,j,k,l;	/* Counter variable */
 int	L=4;	// Minimum length of fragment for FindMaxSub matching
 int	N=15;	// Minimum length of fragment for alignment matching
 int	K1=15;	// Min Stepsize for structural alignment starts 
 int	K2=15;	// Min Stepsize for structural alignment starts 
 int    K=7;    // Max number of steps for structural alignments
 int    Kmax=25;    // Max number of steps for structural alignments
 int    minatoms=19; /* 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          minsim=3.0;
 double          maxpvalue=1,maxpvalue1=1.;
 double          lastscore,dynscore,score,maxrms,maxrms1,maxscore,maxscore1;
 double          pvalue,temp;
 double          bestscore=0.;
 int             maxatoms,maxpos,minpos,atoms,atoms2;
 int             maxatoms1;
 int             maxselect1_0[MAXATMS],maxselect1_1[MAXATMS];
 int             maxselect2_0[MAXATMS],maxselect2_1[MAXATMS];
 int             traceback[MAXATMS];
 int             traceback2[MAXATMS];
 int             besttrace[MAXATMS];
 int             maxsub_15=0;
 int             maxsub_20=0;
 int             maxsub_25=0;
 int             maxsub_30=0;
 int             maxsub_35=0;
 double          bestsizescore[MAXATMS];
 double          a1,b1,identity=0;
 /* Initialize */
 files=0;
 verbose_flag=FALSE;
 output_flag=FALSE;
 out_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],"-pdb")==0)
      {
	out_flag=TRUE;
	i++;
	strcpy(outfile,argv[i]);
      }
    else if (strcmp(argv[i],"-min")==0)
      {
	i++;
	minatoms=atoi(argv[i]);
      }
    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)
       {
	  if (verbose_flag){
	    for (j=0;j<m[0].atoms;j++){ // select some
	      printf("TEST %d %d %d\n",j,m[0].atm[j].resnum,m[1].atm[j].resnum);
	      printf("TEST %d %f %f %f\t%f %f %f\n",j,m[0].atm[j].x,m[0].atm[j].y,m[0].atm[j].z
		     ,m[1].atm[j].x,m[1].atm[j].y,m[1].atm[j].z);
	    }
	   }
	 // this is a test for foo.pdb bar.pdb
	 copy_molecule(&m[1],&m[2]);
	 atoms=m[0].atoms;
	 shortest=min(m[0].atoms,m[1].atoms);
	 
	 if (m[0].atoms/K > K1){ K1=m[0].atoms/K;}
	 if (m[1].atoms/K > K2){K2=m[1].atoms/K;}

	 /*	 if (K1>Kmax){K1=Kmax;}
		 if (K2>Kmax){K2=Kmax;} */

	 /* ************************************************************************
	    First we'll do a structural alignment of the two proteins.
	  ************************************************************************ */
	 for (j=0;j<m[0].atoms;j++){ // select some
	   m[0].atm[j].selected=FALSE;
	   bestsizescore[j]=0;
	   traceback[j]=-1;
	 }
	 for (j=0;j<m[1].atoms;j++){ // select some
	   m[1].atm[j].selected=FALSE;
	 }
	 for (j=0;j<m[0].atoms;j++){ // select some
	   m[0].atm[j].selected=TRUE;
	   //printf("TEST %d %d %d\n",j,m[0].atm[j].resnum,m[1].atm[j].resnum);
	   //	   traceback[j]=m[0].atm[j].resnum-1;
	   //m[1].atm[m[0].atm[j].resnum-1].selected=TRUE;
	   traceback[j]=j;
	   m[1].atm[j].selected=TRUE;
	 }
	 k=0;
	 copy_coord(&m[2],&m[1]);
	 center_molecule(&m[0]);
	 center_molecule(&m[1]);
	 rms=superimpose_molecules(&m[0],&m[1],s,traceback);	  
	 dynscore=Levitt_Gerstein(&m[0],&m[1],traceback2);
	 //printf("dynscore0: %d %d %d   \t%f %f\trmsd:\t%f\n",a,b,k,dynscore,bestscore,rms);
	 //for (i=0;i<m[0].atoms;i++){printf("TRACE-0:  %d  %d\n",i,traceback[i]);}
	 //output_results();

	 //dynscore=dyntrace(&m[0],&m[1],traceback);
	 //rms=superimpose_molecules(&m[0],&m[1],s,traceback);	  
	 //for (i=0;i<m[0].atoms;i++){printf("TRACE-0:  %d  %d\n",i,traceback[i]);}
	 k++;
	 if (dynscore > bestscore){
	   if (verbose_flag) {printf ("Newbestscore1: %f %f\n",dynscore,bestscore);}
	   if (fraction_flag){
	     nalign=0; for (j=0;j<m[0].atoms;j++){if(traceback2[j]>=0){nalign++;}}
	     printf("dynscore0: %d %d %d   \t%f\trmsd:\t%f\t%d\t%f\n",a,b,k,dynscore,rms,nalign,
		  superimpose_molecules(&m[0],&m[1],s,traceback2));
	   }
	   atoms=0;
	   bestscore=dynscore;
	   for (i=0;i<m[0].atoms;i++){
	     besttrace[i]=traceback[i];
	     //printf ("TEST %d %d \n",i,besttrace[i]);
	     if (traceback2[i]>=0){
	       atoms++;
	     }
	   }
	 }
	 //if (FALSE){
	 //printf("dynscore1: %d %d %d   \t%f %f\trmsd:\t%f\n",a,b,k,dynscore,bestscore,rms);
	 for (a=0;a<m[0].atoms-N;a+=K1){
	   for (b=0;b<m[1].atoms-N;b+=K2){
	     //printf ("Aligning:\t%d %d\n",a,b);
  	     for (j=0;j<m[0].atoms;j++){
	       m[0].atm[j].selected=FALSE;
	       traceback[j]=-1;
	     }
	     for (j=0;j<m[1].atoms;j++){
	       m[1].atm[j].selected=FALSE;
	     }
	     for (j=0;j<N;j++){ // select some
	       k=a+j;
	       m[0].atm[k].selected=TRUE;
	       l=b+j;
	       m[1].atm[l].selected=TRUE;
	       traceback[k]=l;
	     }
	     k=0;
	     copy_coord(&m[2],&m[1]);
	     center_molecule(&m[0]);
	     center_molecule(&m[1]);
	     rms=superimpose_molecules(&m[0],&m[1],s,traceback);	  
	     dynscore=dyntrace(&m[0],&m[1],traceback2);
	     //rms=superimpose_molecules(&m[0],&m[1],s,traceback);	  
	     //for (i=0;i<m[0].atoms;i++){printf("TRACE-1:  %d  %d\n",i,traceback[i]);}
	     k++;
	     if (dynscore > bestscore){
	       if (verbose_flag) {printf ("Newbestscore2: %f %f\n",dynscore,bestscore);}
	       if (fraction_flag){
		 nalign=0;for (j=0;j<m[0].atoms;j++){if(traceback2[j]>=0){nalign++;}}
		 printf("dynscore1: %d %d %d   \t%f\trmsd:\t%f\t%d\t%f\n",a,b,k,dynscore,rms,nalign,
		  superimpose_molecules(&m[0],&m[1],s,traceback2));
	       }
	       atoms=0;
	       bestscore=dynscore;
	       for (i=0;i<m[0].atoms;i++){
		 besttrace[i]=traceback2[i];
		 if (traceback2[i]>=0){
		   atoms++;
		 }
	       }
	     }
	     while ( k< 20 && dynscore != lastscore ){
	       for (i=0;i<m[0].atoms;i++){traceback[i]=traceback2[i];}
	       lastscore=dynscore;
	       copy_coord(&m[2],&m[1]);
	       center_molecule(&m[0]);
	       center_molecule(&m[1]);
	       rms=superimpose_molecules(&m[0],&m[1],s,traceback);
	       dynscore=dyntrace(&m[0],&m[1],traceback2);
	       //rms=superimpose_molecules(&m[0],&m[1],s,traceback);
	       //output_results();
	       //for (i=0;i<m[0].atoms;i++){printf("TRACE-2:  %d  %d\n",i,traceback[i]);}
	       //printf("RMS-test2 %f\n",rms);
	       //output_results();
	       if (dynscore > bestscore){
		 if (verbose_flag) {printf ("Newbestscore3: %f %f\n",dynscore,bestscore);}
		 if (fraction_flag){
		   nalign=0;for (j=0;j<m[0].atoms;j++){if(traceback2[j]>=0){nalign++;}}
		   printf("dynscore2: %d %d %d   \t%f\trmsd:\t%f\t%d\t%f\n",a,b,k,dynscore,rms,nalign,
			  superimpose_molecules(&m[0],&m[1],s,traceback2));
		 }
		 atoms=0;
		 bestscore=dynscore;
		 for (i=0;i<m[0].atoms;i++){
		   besttrace[i]=traceback2[i];
		   if (traceback2[i]>=0){
		     atoms++;
		   }
		 }
	       }
	       k++;
	     }
	   }
	   }
	 //{


	 /* ************************************************************************
	    This is for postprocessing of structural alignment 
	    ************************************************************************ */
	 for (j=0;j<m[0].atoms;j++){
	   m[0].atm[j].selected=FALSE;
	 }
	 for (j=0;j<m[1].atoms;j++){
	   m[1].atm[j].selected=FALSE;
	 }
	 for (i=0;i<m[0].atoms;i++){
	   if (besttrace[i]>=0 ){
	     m[0].atm[i].selected=TRUE;
	     m[1].atm[besttrace[i]].selected=TRUE;
	   }
	 }
	 // Identification of fraction similar
	 a=b=1231.23;
	 for (i=0;i<m[0].atoms;i++){
	   a1++;
	   if (strcmp(m[0].atm[i].residue, m[1].atm[besttrace[i]].residue)==0) {b1++;}
	   //	   printf("TEST: %s\t%s\t%f %f %f\n",m[0].atm[i].residue, m[1].atm[besttrace[i]].residue,a1,b1,identity);
	 }
	 identity=b1/a1;
	 if (verbose_flag){
	   for (i=0;i<m[0].atoms;i++){
	     printf ("Struct-align: %d %d \n",i,besttrace[i]);
	   }
	 }
	 copy_coord(&m[2],&m[1]);



	 /* ************************************************************************
	    This is for  testing
	    ************************************************************************ */

	 /*
	 for (i=0;i<m[0].atoms;i++){m[0].atm[i].selected=FALSE;besttrace[i]=-1;}
	 for (i=0;i<m[1].atoms;i++){m[1].atm[i].selected=FALSE;}
	 j=0;
	 for (i=0;i<m[0].atoms;i++){
	   while ( m[0].atm[i].resnum > m[1].atm[j].resnum && j <= m[1].atoms ){j++;}
	   if (m[0].atm[i].resnum ==m[1].atm[j].resnum){
	     m[0].atm[i].selected=TRUE;
	     m[1].atm[j].selected=TRUE;
	     besttrace[i]=j;
	   }
	 }


	 for (i=0;i<m[0].atoms;i++){
	   printf("TRACE:  %d   %d  %d  %d %f\n",i,besttrace[i],m[0].atm[i].resnum,m[1].atm[besttrace[i]].resnum,sqrt(m[0].atm[i].rms));
	 }
	 */
	 /*
	 output_file("out1.pdb");
	 k=0;
	 center_molecule(&m[0]);
	 center_molecule(&m[1]);
	 rms=superimpose_molecules(&m[0],&m[1],s,besttrace);
	 dynscore=dyntrace(&m[0],&m[1],traceback2);
	 //output_file("out2.pdb");
	 //output_results();
	 	 
	 for (i=0;i<m[0].atoms;i++){
	   printf("TRACE:  %d  %d  %d %f\n",i,besttrace[i],traceback2[i],sqrt(m[0].atm[i].rms));
	   if (traceback2[i]>=0){k++;}
	 }
	 a=b=k=0;
	 printf("dynscore: %d %d %d   \t%f\trmsd:\t%f\t%f\n",a,b,k,dynscore,rms,
		  superimpose_molecules(&m[0],&m[1],s,traceback2));
	 
	 exit();
	 */

	 /*
	   for (i=0;i<m[1].atoms;i++){
	   besttrace[i]=i;
	   }
	   atoms=m[1].atoms;
	 */

	 /* ************************************************************************
	    Now it is time to actually find the best subalignment
	    ************************************************************************ */

	 /* ************************************************************************
	    First we'll try to do this the otherway, i.e. adding atoms
	    ************************************************************************ */
	 //printf ("testing %d\n",atoms);

	 maxatoms=0;
	 maxrms=0.0;
	 maxscore=0.0;
	 maxpvalue=1.0;
	 for (i=0;i<atoms-L;i++){
	   for (j=0;j<=m[0].atoms;j++){ // select some
	     m[0].atm[j].selected=FALSE;
	   }
	   for (j=0;j<=m[1].atoms;j++){ // select some
	     m[1].atm[j].selected=FALSE;
	   }
	   k=0;
	   j=i;
	   while(k<L && j<m[0].atoms){ // select some
	     if (besttrace[j]>=0){
	       m[0].atm[j].selected=TRUE;
	       m[1].atm[besttrace[j]].selected=TRUE;
	       //printf("TESTing %d %d\n",j,besttrace[j]);
	       k++;
	     }
	     j++;
	   }
	   j=L;
	   
	   copy_coord(&m[2],&m[1]);
	   center_molecule(&m[0]);
	   center_molecule(&m[1]);
	   //printf("TEST0: %d %f %f\n",i,m[0].atm[0].x,m[1].atm[0].x);
	   rms=superimpose_molecules(&m[0],&m[1],s,besttrace);
	   //printf("TEST3\n");output_results();

	   while (j<atoms && rms< (double)(5*(j+225)/225)){
	     minrmsd=10000.0;
	     minpos=-1;
	     for (k=0;k<m[0].atoms;k++) {
	       //printf("TEST-k: %d\t%f\t%f\t%d\t%d\n",k,minrmsd,m[0].atm[k].rms,m[0].atm[k].selected,besttrace[k]);
	       if (m[0].atm[k].rms < minsim && ! m[0].atm[k].selected ){
		 j++;
		 m[0].atm[k].selected=TRUE;
		 m[1].atm[k].selected=TRUE;
	       }
	       if (m[0].atm[k].rms <= minrmsd && ! m[0].atm[k].selected && besttrace[k]>=0 ){
		 minpos=k;
		 minrmsd=m[0].atm[k].rms;
	       }
	     }
	     j++;
	     //printf("Selecting %d %f %d %d \n",minpos,minrmsd,m[0].atm[minpos].selected,besttrace[minpos]);
	     if(minpos<0){
	       if (verbose_flag){
		 printf ("ERROR1, not selecting any atom: %d\n",j);
		 for (i=0;i<m[0].atoms;i++){
		   printf("SEL-TEST:  %d  %d %d %d\n",i,besttrace[i],m[0].atm[i].selected,m[0].atm[besttrace[i]].selected);
		 }
	       }

	     }
	     else{
	       m[0].atm[minpos].selected=TRUE;
	       m[1].atm[besttrace[minpos]].selected=TRUE;
	     }
	     copy_coord(&m[2],&m[1]);
	     center_molecule(&m[0]);
	     center_molecule(&m[1]);
	     rms=superimpose_molecules(&m[0],&m[1],s,besttrace);	   
	     score=Levitt_Gerstein(&m[0],&m[1]);
	     pvalue=LG_pvalueF(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){
	       //printf  ("TEST\t%d\t%d\%f",j,atoms,fraction);
	       if ((double)j/atoms <= fraction ){
		 printf("FRAC:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
			(double)j/m[0].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 (minatoms < 0){
		 if ((j-shortest) >= minatoms ) {
		   if (fraction_flag){
		     printf("GOOD:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
			    (double)j/m[0].atoms,j,rms,score,
			    pvalue);
		   }
		   //printf ("TEST:\t%d\t%d\t%f\t%f\t%f\n",atoms,j,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;}
		 }
	       }else{
		 if (fraction_flag){
		   printf("GOOD:\t%.2f\t%d\t%.1f\t%.1f\t%e\n",
			  (double)j/m[0].atoms,j,rms,score,
			  pvalue);
		 }
		 //printf ("TEST:\t%d\t%d\t%f\t%f\t%f\n",atoms,j,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("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;


	   if (ma[0]>ma[1]){
	     length=ma[0];
	   }else{
	     length=ma[1];
	   }
           ss=-1000*(log10(maxpvalue)+0.6108)/length;
	   printf("SCORE-1:\t%.2f\t%d\t%.1f\t%.1f\t%e\t%e\n",
		(double)maxatoms/(double)atoms,maxatoms,maxrms,maxscore,maxpvalue,ss);
	   if (align_flag){
	     for (i=0;i<m[0].atoms;i++){m[0].atm[i].selected=maxselect1_0[i]; traceback[i]=-1;}
	     for (i=0;i<m[1].atoms;i++){m[1].atm[i].selected=maxselect1_1[i];}
	     j=0;
	     for (i=0;i<m[0].atoms;i++){
	       if (m[0].atm[j].selected){
		 while(!m[1].atm[j].selected){j++;}
		 traceback[i]=j;
		 j++;
	       }
	     }
	     center_molecule(&m[0]);
	     center_molecule(&m[1]);
	     rms=superimpose_molecules(&m[0],&m[1],s,traceback);	   
	     output_results();
	   }
	   if (out_flag){
	     for (i=0;i<m[0].atoms;i++){m[0].atm[i].selected=maxselect1_0[i]; traceback[i]=-1;}
	     for (i=0;i<m[1].atoms;i++){m[1].atm[i].selected=maxselect1_1[i];}
	     j=0;
	     for (i=0;i<m[0].atoms;i++){
	       //printf("SELE1: %d %d \n",i,m[0].atm[j].selected);
	       if (m[0].atm[j].selected){
		 while(!m[1].atm[j].selected){j++;}
		 //		 printf("SELET: %d %d \n",i,j);
		 traceback[i]=j;
		 j++;
	       }
	     }
	     /*	     for (i=0;i<m[0].atoms;i++){
	       printf("SELE: %d %d %d %d \n",i,m[0].atm[i].selected,traceback[i],m[1].atm[traceback[i]].selected);
	       }	*/ 
	     center_molecule(&m[0]);
	     center_molecule(&m[1]);
	     rms=superimpose_molecules(&m[0],&m[1],s,traceback);	   
	     output_file(outfile);
	   }

	 //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);

	   //printf("IDENTITY:   \t%f\n",identity);
		
	 
	 if (b_flag){
	   for (i=0;i<m[0].atoms;i++){
	     if (bestsizescore[i]>0){
	       printf("SCORE-TEST\t%d\t%f\t%e\n",i,bestsizescore[i],LG_pvalueF(i,bestsizescore[i]));
	     }
	   }
	 }
       }
     else
       printf("Usage: rms [ -a -f -v ] 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,lastresnum;
  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;
	  lastresnum=(int)-9999;
	  //printf("TEst1 %s %d %d\n",name,resnum,lastresnum);
#ifdef  NOZLIB
	  while(fgets(buff,255,fp)!=NULL)
#else
	  while(gzgets(fp,buff,255)!=NULL)
#endif
	    {
	      //printf("TEst2 %s %d %d\n",name,resnum,lastresnum);
	      sscanf(buff,"%s %d %s %s",line_flag,&number,name,residue);
	      //printf("TEst3 %s %d %d\n",name,resnum,lastresnum);
	      if (strcmp("ATOM",line_flag)==0)	/* Is it an ATOM entry? */
		{
		  //printf("TEst4 %s %d %d\n",name,resnum,lastresnum);
		  sscanf(&buff[22],"%d %lf %lf %lf",&resnum,&x,&y,&z);
		  //printf("TEst %s %d %d\n",name,resnum,lastresnum);
		  if (strcmp("CA",name)==0 && resnum > lastresnum)	   /*  Is it an CA atom ? */
		    {
		      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++;
		      lastresnum=resnum;
		    }
		}
	    }
	  m[i].atoms=atoms;
	  ma[i]=atoms;
#ifdef NOZLIB
	  fclose(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(0);
	}
    }
  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
  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++;
	while (!m2->atm[j].selected){j++;}
	if (m1->atm[i].resnum > last1 || m2->atm[j].resnum > last2 ){
	  //printf("Test-gap>\t%d\t%d\t%d\t%d\t%d\n",i,m1->atm[i].resnum,last1,m2->atm[j].resnum,last2);
	  	  numgap++;
	}
	sum+=(1/(1+m1->atm[i].rms/d0));
	//printf("Test>  %d  %d\t%d\t%f\t%f\n",i,j,numgap,sum,m1->atm[i].rms);
	last1=m1->atm[i].resnum;
	last2=m2->atm[j].resnum;
	j++;
      }
    }
  //printf("TEST>\t%d\t%f\t%f\n",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;
  ln=log((double)N);
  if (N<6) {
    expZ=1.;
    return(expZ);
  }
  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> %d\t%f\t%f\t%f\t%f\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,trace)	/* Find RMS superimposition of m1 on m2 */
    molecule 	    *m1,*m2;		/* Molecules to be superimposed */
    double 	    s[3][3];	/* Final transformation matrix */
    int             trace[];
{
    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		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		dist2,dist,x,y,z;		/* Temporary coordinate variables */
    
    //printf("TEST1: %d %f\n",0,m2->atm[0].x);
    
    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 (trace[i]>=0 && m1->atm[i].selected){
	  d1= &(m1->atm[i].x);
	  //printf ("TEST aligning 1b %d %d\n",i,trace[i]);
	  for (j=0;j<3;j++)
	    {
	      d2= &(m2->atm[trace[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;
    natoms=0;
    //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);
	/* 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 ("TEST2: %f %f %f\n",x,y,z);
	  printf ("TEST3: %d %f %f\n",i,error,m1->atm[i].rms); */
      }
    for (i=0;i<m1->atoms;i++){
	x=m1->atm[i].x-m2->atm[trace[i]].x;
	y=m1->atm[i].y-m2->atm[trace[i]].y;
	z=m1->atm[i].z-m2->atm[trace[i]].z;
	m2->atm[trace[i]].rms=m1->atm[i].rms=x*x+y*y+z*z;
	if (trace[i]>=0 && m1->atm[i].selected){
	  error+=m1->atm[i].rms;
	  //printf ("TEST aligning2 %d %d %f\n",i,trace[i],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[trace[i]].x,m2->atm[trace[i]].y,m2->atm[trace[i]].z);
	  printf ("TEST3: %d %f %f\n",i,error,m1->atm[i].rms); */
	  
	}
      }
    //printf("TEST %f\n",error);
    error/=(double)(natoms);
    return(sqrt(error));
}
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);
    //if (verbose_flag) {
	/* First output RMS deviations by atom */
	printf("\n");
	printf("RMS deviations by atom\n");
	printf("i\tAtom\tName\tresidue\tRMS deviation\tj\tAtom\tName\tresidue\tRMS deviation\n");
	j=0;
	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){
	    while(!m[1].atm[j].selected){j++;}
	    printf("%d\t%d\t%s\t%3s %5d\t%8.4e\t%d\t%d\t%s\t%3s %5d\t%8.4e\n",
		 i,m[1].atm[j].number,m[0].atm[i].name,
		 m[0].atm[i].residue,m[0].atm[i].resnum,
		 sqrt(m[0].atm[i].rms),j,
		 m[1].atm[j].number,m[1].atm[j].name,
		 m[1].atm[j].residue,m[1].atm[j].resnum,
		 sqrt(m[1].atm[j].rms));
	    j++;
	  }
	}
	//}
}
output_file(filename)
    char filename[800];
{
    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   %-1d.00     \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,
			m[0].atm[i].selected);
	      }
	    fprintf(fp,"END \n");
	    for (i=0;i<m[0].atoms;i++)
	      {
		fprintf(fp,"ATOM   %-4d  %-4s%-3s  %-4d     %7.3lf %7.3lf %7.3lf  1.0   %-1d.50     \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,
			m[0].atm[i].selected);
	      }
	    fprintf(fp,"END \n");
	}
    fclose(fp);
}

double dyntrace(m1,m2,trace)
    molecule 	*m1,*m2;
    int trace[];
{

  double bestscore,startfact,maxscorex,maxscorey;
  int x, y, type,bestx,besty,i,maxy,maxx;
  double dynscore;
  double skip1, skip2[MAXATMS],minscore;
  double prescore[MAXATMS],  lastscore, bestjump;
  double score, diagtest, test2,tempscore;
  double temp;
  int length,length2,seq[MAXATMS];
  int xminone,yminone;
  int prex[MAXATMS][MAXATMS],prey[MAXATMS][MAXATMS];
  int xlast,ylast[MAXATMS];
  
  length=m1->atoms;
  length2=m2->atoms;
  /*
    c     init dimensions and bestscore
    c     do the bottom row (y=1)
  */
  bestx=0;
  besty=0;
  minscore=0.0;
  bestscore = 0.0;
  skip2[0]=0.0;
  score=prescore[0]=sdist(m1,m2,0,0);
  bestscore = max(bestscore,score);
  bestx=0;
  besty=0;
  prex[0][0]=-1;
  prey[0][0]=-1;
  ylast[0]=0;
  //      C        write (*,*) 'test1:',1,1,prescore(1),' ',prex(1,1),' ',prey(1,1)
  y=0;
  x=0;
  //      printf("TEST2a: %d %d %d %d  %f %f \n",x,y, prex[y][x], prey[y][x],prescore[0],score);
  for (x = 1; x<length;x++){
    prescore[x]=sdist(m1,m2,0,x);
    prex[0][x]=x-1;
    prey[0][x]=-1;
    ylast[x]=0;
    //	printf("TEST2a: %d %d %d %d  %f %f \n",x,y, prex[y][x], prey[y][x],prescore[x],score);
    //C          bestscore = max(bestscore,score)
    if (prescore[x] > bestscore) {
      bestscore=prescore[x];
      bestx=x;
      besty=0;
    }
    skip2[x] = 0.0;
    //C          write (*,*) 'test1:',x,1,prescore(x),' ',prex(1,x),' ',prey(1,x)
  }
  maxscorey=minscore;
  maxy=0;
  skip1=minscore;
  for (y = 1; y<length2;y++){
    yminone=y-1;
    if (maxscorey < score+penopen) {
      maxy=yminone;
      maxscorey=score+penopen;
    }
    //c     for x=0 
    x = 0;
    xlast=0;
    score=sdist(m1,m2,y,x);
    prex[y][x]=-1;
	prey[y][x]=yminone;
	
        lastscore = score;
	//C        bestscore = max(bestscore,score)
        if (score > bestscore) {
          bestscore=score;
          bestx=x;
          besty=y;
	}
        skip1 = minscore;
	//	printf("TEST2b: %d %d %d %d  %f %f \n",x,y, prex[y][x], prey[y][x],sdist(m1,m2,y,x),score);

        for (x = 1; x<length;x++){
          xminone=x-1;
          
          temp = sdist(m1,m2,y,x);
          diagtest = prescore[xminone];
          test2 = skip2[xminone];

	  //	  printf("TEST1:  %d %d %f %f %f %f \n",x,y,diagtest,test2,skip1,temp);
	  //C          score=max(diagtest,max(test2,skip1))+temp
          if (diagtest >= skip1 && diagtest >= test2) {
            score=diagtest+temp;
	    prex[y][x]=xminone;
	    prey[y][x]=yminone;
	  }else if (test2 >= skip1) {
	    score=test2+temp;
            prex[y][x]=xminone;
            prey[y][x]=ylast[xminone];
          }else{
            score=skip1+temp;
            prex[y][x]=xlast;
            prey[y][x]=yminone;
	  }
	  //C          bestscore = max(bestscore,score)
          if (score > bestscore) {
            bestscore=score;
            bestx=x;
            besty=y;
	  }

          bestjump = prescore[xminone]+penopen;

          if (bestjump > skip1) {
            skip1 =bestjump;
            xlast=xminone;
	  }
	  //C          skip2(x) = max(bestjump,skip2(xminone))
          if (skip2[xminone] < bestjump) {
            skip2[xminone] = bestjump;
            ylast[xminone]=yminone;
	  }

          prescore[xminone] = lastscore;
          lastscore = score;
	  //C          write (*,*) 'test3:',x,y,prescore(xminone),score,' ',prex(y,x),' ',prey(y,x)
	  //	  printf("TEST2c:  %d %d %d %d  %f %f \n",x,y, prex[y][x], prey[y][x],temp,score);
	}
        prescore[0]=sdist(m1,m2,y,0);
      }

      if (maxscorey < score) {
        maxy=length2-1;
        maxscorey=score;
      }

      maxscorex=score;
      maxx=length-1;
      for (x = length-1;x>=0;x--){
        if (maxscorex < prescore[x]) {
          maxx=x;
          maxscorex=prescore[x];
	}
      }

      //            printf("TRACE-test> %f %f %f %d %d\n",score,maxscorey,maxscorex,maxx,maxy);
      dynscore=bestscore;
      
      //      dynscore=max(score,max(maxscorey,maxscorex));
      
	/*      if (score > maxscorey && score > maxscorex) {
		dynscore=score;
		bestx=length-1;
		besty=length2-1;
		}else if ( maxscorey > maxscorex) {
		dynscore=maxscorey;
		bestx=length-1;
		besty=maxy;
		}else{
		dynscore=maxscorex;
		bestx=maxx;
		besty=length2-1;
		}
	*/

      for (x=0;x<length;x++){
        trace[x]=-1;
      }

      y=besty;
      x=bestx;
      bestscore=0;
      //printf("TEST-78 %d %d\n",besty,bestx);
      while (y > -1 && x > -1){
	//printf("TRACE-ds: %d %d %d %d %f\n",x,y,prex[y][x],prey[y][x],sdist(m1,m2,y,x));
	//bestscore+=sdist(m1,m2,y,x);
	//if (x-prex[y][x]>1 ||y-prey[y][x]>1 ){bestscore+=penopen;}

	trace[x]=y;
	i=x;
	x=prex[y][x];
	y=prey[y][i];
      }
      for (x=0;x<length;x++){m1->atm[x].selected=FALSE;}
      for (x=0;x<length2;x++){m2->atm[x].selected=FALSE;}
      for (x=0;x<m1->atoms;x++){
	if(trace[x]>=0){
	  m1->atm[x].selected=TRUE;
	  m2->atm[trace[x]].selected=TRUE;
	}
      }
      //printf ("Align-test %f %f %f\n",dynscore,bestscore,dynscore-bestscore);
      //for (x=0;x<m1->atoms;x++){printf("TRACE-ali:  %d  %d\n",x,trace[x]);}
      return(dynscore);
}


double sdist(m1,m2,j,i)
    molecule 	*m1,*m2;
    int i,j;
{
  double score;
  double dist;
  dist=
    (m1->atm[i].x-m2->atm[j].x)*(m1->atm[i].x-m2->atm[j].x)+
    (m1->atm[i].y-m2->atm[j].y)*(m1->atm[i].y-m2->atm[j].y)+
    (m1->atm[i].z-m2->atm[j].z)*(m1->atm[i].z-m2->atm[j].z);
  score=M/(1+dist/d0);
  /*
    printf("SDIST:\t%d\t%d\t%f\t%f\n",i,j,sqrt(dist),score);
    printf("XYZ-1:\t%f\t%f\t%f\n",m1->atm[i].x,m1->atm[i].y,m1->atm[i].z);
    printf("XYZ-2:\t%f\t%f\t%f\n",m2->atm[j].x,m2->atm[j].y,m2->atm[j].z);
  */
  return(score);
}

/*double max(i,j)
     double i,j;
{
  if (i>j){
    return(i);
  }else{
    return(j);
  }
}

double min(i,j)
     double i,j;
{
  if (i>j){
    return(j);
  }else{
    return(i);
  }
}
*/

copy_coord(m1,m2)
     molecule 	*m1,*m2;
{
  int i;
  for (i=0;i<m1->atoms;i++){
    m2->atm[i].x=m1->atm[i].x;  
    m2->atm[i].y=m1->atm[i].y;  
    m2->atm[i].z=m1->atm[i].z;  
  }
}

copy_molecule(m1,m2)
     molecule 	*m1,*m2;
{
  int i;
  m2->atoms=m1->atoms;
  m2->xcen=m1->xcen;
  m2->ycen=m1->ycen;
  m2->zcen=m1->zcen;
  strcpy(m1->filename,m2->filename);
  for (i=0;i<m1->atoms;i++){
    copy_res(m1,m2,i,i);
  }
}

copy_res(m1,m2,i,j)
    molecule 	*m1,*m2;
    int i,j;
{
  m2->atm[j].x=m1->atm[i].x;  
  m2->atm[j].y=m1->atm[i].y;  
  m2->atm[j].z=m1->atm[i].z;  
  m2->atm[j].rms=m1->atm[i].rms;  
  //strcpy(m1->atm[i].name,m2->atm[j].name);
  //strcpy(m1->atm[i].residue,m2->atm[j].residue);
  m2->atm[j].number=m1->atm[i].number;  
  m2->atm[j].resnum=m1->atm[i].resnum;
  m2->atm[j].selected=m1->atm[i].selected;  
}


////////////////////////
//The following is added by Fang
double LG_pvalueF(N,score)
     int N;
     double score;
{
    double Z,s1,mean,SD,sc,expZ;
    s1=(double)N;
    mean=pow(s1,0.3264)/(0.0437*pow(s1,0.0003)+0.0790);
    SD=s1/(0.0417*s1+3.3700);
    Z=(score/10-mean)/SD;
    expZ=exp((double)-1.*Z);
    if (Z>20)s1=expZ;
    else s1=1-exp(-expZ);
    return s1;
}
//End
///////////////////////////