/* combine.c++  Routines for combining 2 or more runs of maid.
Procedure:
   Input a set of fit files: stored in file combinefile. Example:
2
finalfit_fuse
fitb.29jun00

   1) Reads all the files, store fit in usual fit[i=1,,,,] and
        fitnum[0] = total number of fits.
   2) Writes this to standard maid.fitfile - as usual.
THEN  Detetes and modifies fits in series of steps:
  3) Looks at all unassigned fits -
       1)Removes identical fits - if one completely overlaps another.
          (This works fine).
       2)If do not completely overlap - want to extend one so it does.
           THIS HAS NOT BEEN WRITTEN YET BECAUSE NO EXAMPLE TO DEBUG.
           Then - takes the longer fit - adds any missing residues from
             shorter fit, starting with known phi/psi and refining.
            The shorter fit is then deleted.
   4)
   5) At this point - have assigned residue numbers to all possible fits.
   6) Next - run combinefit:
       This looks through fits to find fits with some overlap of seq. assighment.
*/
#include "drawwindow.h"
#include "math.h"
#include "trace.h" // includes a number of "extern .... " lines
#include "peptide.h"
#include "pdbtolev.h"
#include "assignseq.h"
#include "symmate.h"
#include "combine.h"
#include "extendfit.h"
#include "sphere.h"
#include "fittopdb.h"
#include "string.h"
#include "stdlib.h" 

extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern struct maiddatafile maidfile;
extern struct spheredata sp[NSPHERES];
extern int fitnum[NFIT];// see trace.c++ for definition
extern int verbose;
extern int side[AA][NBRN][NPTS];  // defined globally (originally in pdbtolev.c++)  because needed in assignseg.c++ 

#define NEND 50
FILE *fperrmes; // file where write error message

void readonefitfile(struct geometry& geo,char fitfile[],int istart)
/* Similar to readfit except that fits are put in fit[istart,istart+1,..] 
*/
{
   int i,irank0,nfits;
   FILE *fp;

   if( (fp=fopen(fitfile,"r")) == NULL){
       fitnum[0] = 0; // indicates no fit residues
       printf("ERROR: fitfile %s does not exist\n",fitfile);
       exit(0);
   }
   else{ //fitfile opened
       fscanf(fp,"%d",&nfits); // total  # of fits;
       if(verbose>=1)
          printf("Reading %d fits from fitfile %s\n",nfits,fitfile);
       for(i=1;i<=nfits;i++){ // each fit
          irank0=0;  // there are no rankings for files read in readfit
          readonefit(fp,fit,irank0,istart+i-1);
       }// each fit 
       fitnum[0] = fitnum[0]+nfits;
       fclose(fp);  
   } // fitfile exists 
   geo.fuse[0] = geo.fuse[1] = geo.fuse[2] = -1;  //Default values
}

void readcombinefile(struct geometry& geo,char combinefile[])
{
   int i,ii,icomb,istart;
   FILE *fpin;
   char fitfile[300];

   if( (fpin=fopen(combinefile,"r")) == NULL){
     printf("ERROR: combinefile does not exist\n" );
     exit(0);
   }
   else{ //combinefile opened
     fscanf(fpin,"%d",&icomb); // total  # of files;
     if(verbose>=1)
        printf("Reading %d fitfiles from file %s\n",icomb,combinefile);
     fitnum[0]=0;
     for(i=1;i<=icomb;i++){ // read eachfile
        fscanf(fpin,"%s",fitfile); // total  # of files;
        istart = fitnum[0]+1;
        readonefitfile(geo,fitfile,istart); // assign firstfit to fitnum[0]+1 = istart
        for(ii=istart;ii<=fitnum[0];ii++){// set fit.val2 = +/- icomb
           if(fit[ii].val2<0)
              fit[ii].val2 = -i;
           else
              fit[ii].val2 = i;    
        }        
     }
     savefit(); // write the new combined fit file
   }// combinefile read
   fclose(fpin);  
}

void analyzlap(int ifit,int ifit2,int ires,int res[],float lap[10][10],
      int *laptype,int *lapres)
/*
   Analyses all overlaps of residue ires of fit[ifit] with fit ifit2
      res[1,2,..] are the residue of ifit2 that overlap ires
      res[0] = number of overlaps
      lap[1,2,..][1,2,4] is overlap dist of Ca, C and N of residues that overlap
   checks the over lap for one residue - decides if they are identical,
     or just meet at end, or indicate some problem
  laptype = 1: standard identical overlap
  laptype = 2: cterm to nterm
  laptype = 3: nterm to nterm  BAD
  laptype = 4: nterm to cterm
  laptype = 5: cterm to cterm  Bad
  laptype = 10: something funny
   lap[i][1,2,4] = overlap dist for Ca,C and N
*/
{
   int i;
   float distc,distn;

   *laptype = 0; // default 
   for(i=1;i<=res[0];i++){ // each residue that overlaps
      if(lap[i][1]< 1.0){ //Close Ca
        if(ires==0){ // nterm residue
          distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
          if(distc<1.0){
            if(*laptype ==1)
               printf("Error, 2 identical residues in anlayselap\n");
            else{
               *laptype = 1; //standard identical overlap
               *lapres = res[i];
            }
          }
          else if(res[i]==fit[ifit2].nres-1)
               *laptype = 2; // cterm to nterm meet
          else if(res[i]==0)
               *laptype = 3; // nterm to nterm meet - BAD
          else{
               *laptype = 10;
               *lapres = res[i];
          } 
        } // ires = 0
        else if(res[i]==0){ // nterm residue
          distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
          if(distc<1.0){
            if(*laptype ==1)
               printf("Error, 2 identical residues in anlayselap\n");
            else{
               *laptype = 1; //standard identical overlap
               *lapres = res[i];
            }
          }
          else if(ires==fit[ifit].nres-1)
               *laptype = 2; // cterm to nterm meet
          else if(ires==0)
               *laptype = 3; // nterm to nterm meet - BAD
          else{
               *laptype = 10; 
               *lapres = res[i];
          } 
        } // ires = 0
        else if(ires == fit[ifit].nres-1){ // cterm
            distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
            if(distn<1.0 ){
                if(*laptype ==1)
                   printf("Error, 2 identical residues in anlayselap\n");
                else{
                   *laptype = 1; //standard identical overlap
                   *lapres = res[i];
                }
            }
            else if(res[i]==fit[ifit2].nres-1)
               *laptype = 5; // cterm to cterm meet BAD
            else if(res[i]==0)
               *laptype = 4; // nterm to cterm meet 
            else{
               *laptype = 10; 
               *lapres = res[i];
            } 
        } // ires = 0 
        else if(res[i] == fit[ifit2].nres-1){ // cterm
            distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
            if(distn<1.0 ){
                if(*laptype ==1)
                   printf("Error, 2 identical residues in anlayselap\n");
                else{
                   *laptype = 1; //standard identical overlap
                   *lapres = res[i];
                }
            }
            else if(ires==fit[ifit].nres-1)
               *laptype = 5; // cterm to cterm meet BAD
            else if(ires==0)
               *laptype = 4; // nterm to cterm meet 
            else
               *laptype = 10; 
         } // ires = 0 
         else{ // central res    
           distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
           distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
           if( (distn<1.0)&&(distc<1.0) ){
              if(*laptype ==1)
                 printf("Error, 2 identical residues in anlayselap\n");
              else{
                *laptype = 1; //standard identical overlap
                *lapres = res[i];
              }
           }// distn and distc OK
           else
               *laptype = 10; 
        } // central res
      } // close Ca
   } // i                
}


void analyzbad(int *iquit,int ifit,int ifit2,int ires,int res[],float lap[10][10],
      int ibad[MAXIBAD][IBADRES][2])
/*
   Finds bad contacts - 
     ASSUMES HAVE ALLREADY ELIMINATED IDENTICAL FITS - THUS ANY OTHER CONTACTS
       ARE EITHER MEETING OF ENDS - OR ARE BAD
   Analyses all overlaps of residue ires of fit[ifit] with fit ifit2
      res[1,2,..] are the residue of ifit2 that overlap ires
      res[0] = number of overlaps
      lap[1,2,..][1,2,4] is overlap dist of Ca, C and N of residues that overlap
   checks the over lap for one residue - decides if they are identical,
     or just meet at end, or indicate some problem
  laptype = 1: standard identical overlap
  laptype = 2: cterm to nterm
  laptype = 3: nterm to nterm  BAD
  laptype = 4: nterm to cterm
  laptype = 5: cterm to cterm  Bad
  laptype = 10: something funny
   lap[i][1,2,4] = overlap dist for Ca,C and N
*/
{
   int i,ii,laptype;
   float distc,distn;

   *iquit=0;  // default OK
   laptype = 10; // default  - assume bad
   for(i=1;i<=res[0];i++){ // each residue that overlaps
      if(lap[i][1]< 1.0){ //Close Ca - look for ends meeting
        if(ires==0){ // nterm residue
          distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
          if(distc>1.0){
             if(res[i]==fit[ifit2].nres-1){
                 laptype = 2; // cterm to nterm meet
                 if(verbose>=1)
                    printf(" nterm of Fit %d  meets cterm of fit %d \n",ifit,ifit2);
             }
             else if(res[i]==0)
                 laptype = 3; // nterm to nterm meet - BAD
          }//distc>1
        } // ires = 0
        else if(res[i]==0){ // nterm residue
          distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
          if(distc>1.0){
              if(ires==fit[ifit].nres-1){
                  laptype = 2; // cterm to nterm meet
                  if(verbose>=1)
                     printf(" nterm of Fit %d  meets cterm of fit %d \n",ifit,ifit2);
             }
               else if(ires==0)
                   laptype = 3; // nterm to nterm meet - BAD
          }// distc
        } // ires = 0
        else if(ires == fit[ifit].nres-1){ // cterm
            distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
            if(distn>1.0 ){
               if(res[i]==fit[ifit2].nres-1)
                  laptype = 5; // cterm to cterm meet BAD
               if(res[i]==0){
                  laptype = 4; // nterm to cterm meet 
                  if(verbose>=1)
                     printf(" Cterm of Fit %d meets nterm of fit %d \n",ifit,ifit2);
               }
            }
        } // ires = 0 
        else if(res[i] == fit[ifit2].nres-1){ // cterm
            distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
            if(distn>1.0 ){
                if(ires==fit[ifit].nres-1)
                   laptype = 5; // cterm to cterm meet BAD
                if(ires==0){
                   laptype = 4; // nterm to cterm meet 
                   if(verbose>=1)
                     printf(" Cterm of Fit %d meets nterm of fit %d \n",ifit,ifit2);
                }
                else if(laptype==4)
                   if(verbose>=1)
                      printf(" Cterm of Fit %d meets nterm of fit %d \n",ifit,ifit2);
             } 
         } // ires = 0 
      } // close Ca
      if( (laptype==5)|(laptype==3)||(laptype==10) ){ // if bad
         for(ii=1;ii<=ibad[0][0][0];ii++){
            if(ibad[ii][0][0]==ifit2)// previous bad res from this fit - add to llist
             goto addi;
         }
         ii=ibad[0][0][0] = ibad[0][0][0]+1; // if reach here  -first residue from this fit that is identical
         if(ii>= MAXIBAD-1){
               printf("ERROR - Too many bad fits - quitting combineonefit\n");
               *iquit = 1;
               goto quitrout;
         }
         ibad[ii][0][0] = ifit2;
         addi:;
         ibad[ii][0][1] = ibad[ii][0][1]+1;
         if(ibad[ii][0][1] >=IBADRES-1){
                printf(" ERROR- Too many bad resudue - quitting combineonefit\n");
                *iquit=1;
                goto quitrout;
         }
         ibad[ii][ ibad[ii][0][1] ][0] =ires; 
         ibad[ii][ ibad[ii][0][1] ][1] =res[i];
      }// if bad
   } // i                
   quitrout:;
}


void analyzident(int *iquit,int ifit,int ifit2,int ires,int res[],float lap[10][10],
      int ident[IDENTFITS][MAXIDENT][2])
/*
NOTE - 13feb01 - this return only laptype = 0 or 1 
 I do not know why I did not use every thing in analyzelap

   Analyses all overlaps of residue ires of fit[ifit] with fit ifit2
      res[1,2,..] are the residue of ifit2 that overlap ires
      res[0] = number of overlaps
      lap[1,2,..][1,2,4] is overlap dist of Ca, C and N of residues that overlap
   checks the over lap for one residue - decides if they are identical,
     or just meet at end, or indicate some problem
  laptype = 1: standard identical overlap
  laptype = 2: cterm to nterm
  laptype = 3: nterm to nterm  BAD
  laptype = 4: nterm to cterm
  laptype = 5: cterm to cterm  Bad
  laptype = 10: something funny
   lap[i][1,2,4] = overlap dist for Ca,C and N
*/
{
   int i,laptype,lapres;
   float distc,distn;

   *iquit = 0; // default - OK
   laptype = 0; // default 
   for(i=1;i<=res[0];i++){ // each residue that overlaps
      if(lap[i][1]< 1.0){ //Close Ca
        if(ires==0){ // nterm residue
          distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
          if(distc<1.0){
            if(laptype ==1)
               printf("Error, 2 identical residues in anlayselap\n");
            else{
               laptype = 1; //standard identical overlap
               lapres = res[i];
            }
          }
        } // ires = 0
        else if(res[i]==0){ // nterm residue
          distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
          if(distc<1.0){
            if(laptype ==1)
               printf("Error, 2 identical residues in anlayselap\n");
            else{
               laptype = 1; //standard identical overlap
               lapres = res[i];
            }
          }
        } // ires = 0
        else if(ires == fit[ifit].nres-1){ // cterm
            distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
            if(distn<1.0 ){
                if(laptype ==1)
                   printf("Error, 2 identical residues in anlayselap\n");
                else{
                   laptype = 1; //standard identical overlap
                   lapres = res[i];
                }
            }
        } // ires = 0 
        else if(res[i] == fit[ifit2].nres-1){ // cterm
            distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
            if(distn<1.0 ){
                if(laptype ==1)
                   printf("Error, 2 identical residues in anlayselap\n");
                else{
                   laptype = 1; //standard identical overlap
                   lapres = res[i];
                }
            }
         } // ires = 0 
         else{ // central res    
           distc = distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist
           distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);// N-N dist
           if( (distn<1.0)&&(distc<1.0) ){
              if(laptype ==1)
                 printf("Error, 2 identical residues in anlayselap\n");
              else{
                laptype = 1; //standard identical overlap
                lapres = res[i];
              }
           }// distn and distc OK
        } // central res
      } // close Ca
   } // i                
   if(laptype==1){
      for(i=1;i<=ident[0][0][0];i++){
          if(ident[i][0][0]==ifit2)// previous identical res from this fit - add to llist
             goto addi;
      }
      i=ident[0][0][0] = ident[0][0][0]+1; // if reach here  -first residue from this fit that is identical
      if(i>=IDENTFITS-1){
               printf("ERROR at 1- Too many identical fits - quitting analyziden\n");
               *iquit = 1;
               goto quitrout;
      }
      ident[i][0][0] = ifit2;
      addi:;
      ident[i][0][1] = ident[i][0][1]+1;
//printf("ident[i][0][1] =%d,ires=%d,lapres=%d\n",ident[i][0][1],ident[i][ ident[i][0][1] ][0],ident[i][ ident[i][0][1] ][1]);
      if(ident[i][0][1] >=MAXIDENT-1){
                printf(" ERROR  at 2-  Too many identical residues- quitting analyzeident\n");
               *iquit = 1;
                goto quitrout;
      }
      ident[i][ ident[i][0][1] ][0] =ires; 
      ident[i][ ident[i][0][1] ][1] =lapres;
   }// laptype
   quitrout:;
}


void analyzident2(int *iquit,int ifit,int ifit2,int ires,int res[],float lap[10][10],
      int ident[IDENTFITS][MAXIDENT][2],float lapdist)
/*
NOTE - 13feb01 - Only tests for overlap
*/
{
   int i,laptype,lapres;
   float distc,distn;

   *iquit = 0; // default - OK
   laptype = 0; // default 
   for(i=1;i<=res[0];i++){ // each residue that overlaps
      if(ires>0){ // Necessary for overlap
         distc=  distv(fit[ifit].pos[ires][2],fit[ifit2].pos[res[i]][2]);// C-C dist  
         distn = distv(fit[ifit].pos[ires][4],fit[ifit2].pos[res[i]][4]);//  N-N dist
         if( (distc<lapdist)&&(distn<lapdist) ){// Then overlap
               laptype = 1; //standard identical overlap
               lapres = res[i];
         }
      } // ires>0
   } // i                
   if(laptype==1){
      for(i=1;i<=ident[0][0][0];i++){
          if(ident[i][0][0]==ifit2)// previous identical res from this fit - add to llist
             goto addi;
      }
      i=ident[0][0][0] = ident[0][0][0]+1; // if reach here  -first residue from this fit that is identical
      if(i>=IDENTFITS-1){
               printf("ERROR at 1 in ident2 - Too many identical fits - quitting analyziden\n");
               *iquit = 1;
               goto quitrout;
      }
      ident[i][0][0] = ifit2;
      addi:;
      ident[i][0][1] = ident[i][0][1]+1;
//printf("ident[i][0][1] =%d,ires=%d,lapres=%d\n",ident[i][0][1],ident[i][ ident[i][0][1] ][0],ident[i][ ident[i][0][1] ][1]);
      if(ident[i][0][1] >=MAXIDENT-1){
                printf(" ERROR at 2 in ident2-  Too many identical residues- quitting analyzeident\n");
               *iquit = 1;
                goto quitrout;
      }
      ident[i][ ident[i][0][1] ][0] =ires; 
      ident[i][ ident[i][0][1] ][1] =lapres;
   }// laptype
   quitrout:;
}

void overlap(int ifit,int ifit2,float lapdist,int ires,
         int res[],float lap[10][10])
/* Returns res[0] residues that overlap at dist < lapdist
    lap[i][j] = overlap dist j = 1,2,4 (Ca,C, N atoms) for ith overlap.
*/
{
    int i,i2,nlap,ilap,ii;
    float dist;
         
          nlap = 0;  
          for(i2=0;i2<fit[ifit2].nres;i2++){//Use terminal atoms
             ilap = 0;
             for(i=1;i<=4;i++)
                lap[nlap+1][i]= 10.0; //defautlt indicating non-overlap
             for(ii=1;ii<=4;ii++){// mainchain atoms
               if(ii!=3){ // Not O
                 if((ii==4)&&(i2==0) )
                   goto nextii;
                 if( (ii==2)&&(i2==fit[ifit2].nres-1)  )
                   goto nextii;
                 dist = distv(fit[ifit].pos[ires][1],fit[ifit2].pos[i2][ii]);                 
                 if(dist<lapdist){
                     if(ilap==0){
                        nlap = nlap+1;
                        if(nlap> 9){
                           printf("ERROR: nlap=%d is too big in combineonefit\n",nlap);
                           exit(0);
                        }
                     }
                     ilap = 1;
                     res[nlap] = i2;
                     lap[nlap][ii] = dist;
#if 0
                    printf("fit %d res %d atom %d overlaps fit %d res %d (dist = %5.3f)\n",
                       ifit2,i2,ii,ifit,ires,dist);
#endif
                 }
               } // ii !=3 mainchain
               nextii:;
             }//ii
           }// i2
           res[0]=nlap;
}

void findidentical(int *iquit,int del[],int iass,int ifit,int ident[IDENTFITS][MAXIDENT][2])
/* For each Ca in each fit, test all other fits to see if mainchain
   atoms overlap  (within 1.0 angstroms)
   iass = 0  indicates testing only unassigned fits versus all - and only fit[j>ifit]
   iass = 1  indicates testing only assigned fits versus assigned  - and only fit[j>ifit]
   iass = -1  check all
   Then tests this overlap to see if identical
Returns ident[i][j]:
   ident[0][0][0] = number of different fits ifit has identical residues
   ident[1][0][0] = identical fit 1
   ident[1][0][1] = number of residues in fit 1 that are identical
   ident[1][i][0,1] = ith identical residue in fit 1 corresponds to  residue [0] in ifit and residue [1] in fit 1
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ii,ires,j;
   int res[10];
   float lap[10][10],lapdist;

     ident[0][0][0]=0;
     for(i=0;i<IDENTFITS;i++)
        ident[i][0][1] =0;
     for(ires=0;ires<fit[ifit].nres;ires++){ // Each Ca in fit[i]
       for(j=1;j<=fitnum[0];j++){
          for(ii=1;ii<=del[0];ii++)
             if(del[ii]==j)
               goto nextj;
          if((iass==0)&&(fit[j].val1>=0)&&(j<=ifit) ) // This prevents duplicatiion
              goto nextj;
          if(iass==1){//indicates testing only assigned fits versus assigned
             if( (fit[j].val1>=0)||(j<=ifit) )
                goto nextj;
          }
          lapdist = 1.0;
          overlap(ifit,j,lapdist,ires,res,lap);
          if(res[0]>0){
                analyzident(iquit,ifit,j,ires,res,lap,ident);
                if(*iquit==1)
                   goto quitrout;
          }//res[0] >0
          nextj:;
       }//j
   } // ires
   quitrout:;
}

void findoneidentical(int *iquit,int del[],int iass,int ifit,int ident[IDENTFITS][MAXIDENT][2])
/*
if iass = 1 Only test assigned fits.
if iass = 0  test all fits.
if iass =2  Only test unassigned fits
Same as findidentical - except applied to one fit at a time so try all fits
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ii,ires,j;
   int res[10];
   float lap[10][10],lapdist;

     ident[0][0][0]=0;
     for(i=0;i<IDENTFITS;i++)
        ident[i][0][1] =0;
     for(ires=0;ires<fit[ifit].nres;ires++){ // Each Ca in fit[i]
       for(j=1;j<=fitnum[0];j++){
//  for(j=77;j<=77;j++){
          for(ii=1;ii<=del[0];ii++)
             if(del[ii]==j)
               goto nextj;
          if(j==ifit)
              goto nextj;
          else if( (iass==1)&&(fit[j].val1>=0) )//indicates testing  fits versus only assigned
              goto nextj;
          else if( (iass==2)&&(fit[j].val1<0) )//indicates testing  fits versus only unassigned
              goto nextj;
          lapdist = 1.3;
          overlap(ifit,j,lapdist,ires,res,lap); // test all atoms in fit[j] for overlap with ires in fit[ifit]
          if(res[0]>0){
                analyzident(iquit,ifit,j,ires,res,lap,ident);
                if(*iquit==1)
                   goto quitrout;
          }//res[0] >0
          nextj:;
       }//j
   } // ires
   quitrout:;
}

#if 1
void findoneidentical2(int *iquit,int del[],int iass,int ifit,int ident[IDENTFITS][MAXIDENT][2])
/*
Simple version of findoneidentical:
  Call analyzident2;
   Allows slightly larger lapdist = 1.2
if iass = 1 Only test assigned fits.
if iass = 0  test all fits.
if iass =2  Only test unassigned fits
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ii,ires,j;
   int res[10];
   float lap[10][10],lapdist;

     ident[0][0][0]=0;
     for(i=0;i<IDENTFITS;i++)
        ident[i][0][1] =0;
     for(ires=0;ires<fit[ifit].nres;ires++){ // Each Ca in fit[i]
       for(j=1;j<=fitnum[0];j++){
//  for(j=77;j<=77;j++){
          for(ii=1;ii<=del[0];ii++)
             if(del[ii]==j)
               goto nextj;
          if(j==ifit)
              goto nextj;
          else if( (iass==1)&&(fit[j].val1>=0) )//indicates testing  fits versus only assigned
              goto nextj;
          else if( (iass==2)&&(fit[j].val1<0) )//indicates testing  fits versus only unassigned
              goto nextj;
          lapdist = 1.2;
          overlap(ifit,j,lapdist,ires,res,lap); // tries all atoms in fit[j] return ires[1] = residue with overlap
          if(res[0]>0){
                analyzident2(iquit,ifit,j,ires,res,lap,ident,lapdist);
                if(*iquit==1)
                   goto quitrout;
          }//res[0] >0
          nextj:;
       }//j
   } // ires
   quitrout:;
}
#endif

void findbad(int *iquit,int iass,int ifit,int ibad[MAXIBAD][IBADRES][2])
/* For each Ca in each fit, test all other fits to see if mainchain
   atoms overlap  (within 3.0 angstroms)
   Then tests this overlap to see if bad (have already removed identical fits
   iass = 0  indicates testing only unassigned fits versus all
Returns ident[i][j]:
   ibad[0][0][0] = number of different fits ifit have residues
   ibad[1][0][0] = bad fit 1
   ibad[1][0][1] = number of residues in fit 1 that are bad
   ibad[1][i][0,1] = ith bad residue in fit 1 corresponds to  residue [0] in ifit and residue [1] in fit 1
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ires,j;
   int res[10];
   float lap[10][10],lapdist;

     ibad[0][0][0]=0;
     for(i=0;i<20;i++)
        ibad[i][0][1] =0;
     for(ires=0;ires<fit[ifit].nres;ires++){ // Each Ca in fit[i]
       for(j=1;j<=fitnum[0];j++){
          if((iass==0)&&(fit[j].val1>=0)&&(j<=ifit) ) // This prevents duplicatiion
            goto nextres;
          lapdist = 3.0;
          overlap(ifit,j,lapdist,ires,res,lap);
          if(res[0]>0){
                analyzbad(iquit,ifit,j,ires,res,lap,ibad);
          }//res[0] >0
          nextres:;
       }//j
       if(*iquit==1)
          goto quitrout;
   } // ires
   quitrout:;
}

void combineonefit(int ifit,int ident[IDENTFITS][MAXIDENT][2],int badcont[10])
/* 
NOTE - FOR NOW THIS IS NOT USED
For each Ca in each fit, test all other fits to see if mainchain
   atoms overlap  (within 3.0 angstroms)
Returns ident[i][j]:
   ident[0][0][0] = number of different fits ifit has identical residues
   ident[1][0][0] = identical fit 1
   ident[1][0][1] = number of residues in fit 1 that are identical
   ident[1][i][0,1] = ith identical residue in fit 1 corresponds to  residue [0] in ifit and residue [1] in fit 1
Returns badcont[i]:
   badcont[0][0][0] = number of different fits ifit has bad contact residues
   badcont[1][0][0] = bad contact fit 1
   badcont[1][0][1] = number of residues in fit 1 that are bad contact
   badcont[1][i][0,1] = ith bad contact residue in fit 1 corresponds to  residue [0] in ifit and residue [1] in fit 1
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ires,j;
   int res[10],laptype,lapres;
   float lap[10][10],lapdist;

     ident[0][0][0]=badcont[0]=0;
     for(i=0;i<IDENTFITS;i++)
        ident[i][0][1] =0;
     for(ires=0;ires<fit[ifit].nres;ires++){ // Each Ca in fit[i]
       for(j=1;j<=fitnum[0];j++){
         if(j!=ifit){
           lapdist = 3.0;
           overlap(ifit,j,lapdist,ires,res,lap);
           if(res[0]>0){
                analyzlap(ifit,j,ires,res,lap,&laptype,&lapres);
                if(laptype==1){
//                    printf("Fit %d res %d is idential to fit %d res %d\n",ifit,ires,j,lapres);
                    for(i=1;i<=ident[0][0][0];i++){
                      if(ident[i][0][0]==j)// previous identical res from this fit - add to llist
                         goto addi;
                    }
                    i=ident[0][0][0] = ident[0][0][0]+1; // if reach here  -first residue from this fit that is identical
                    ident[i][0][0] = j;
                    if(i>= IDENTFITS-1){
                      printf("ERROR - ERROR- Too many identical fits  - quitting combineonefit\n");
                      goto quitrout;
                    }
                    addi:;
                    ident[i][0][1] = ident[i][0][1]+1;
                    if(ident[i][0][1] >=MAXIDENT-1){
                        printf(" ERROR- Too many identical resudue - quitting combineonefit\n");
                        goto quitrout;
                    }
                    ident[i][ ident[i][0][1] ][0] =ires; 
                    ident[i][ ident[i][0][1] ][1] =lapres;
                }
                else if( (laptype ==3)||(laptype==5)||(laptype==10) ){
#if 1
                   for(i=1;i<=badcont[0];i++){
                      if(badcont[i]==j) // already noted bad contact - do nothing
                         goto nextres;
                   }
#endif
#if 0
                   for(i=1;i<=badcont[0][0][0];i++){
                      if(badcont[i][0][0]==j) // already noted bad contact - do nothing
                         badcont[i][0][1] = badcont[i][0][1]+1;
                         badcont[i][ badcont[i][0][1] ][0]=ires;
                         badcont[i][ badcont[i][0][1] ][0]=lapres;
                         goto nextres;
                   }
#endif
                   badcont[0]=badcont[0]+1; // if reach here - first bad contact with this fit
                   if(badcont[0]>=10){
                        printf("ERROR: too many bad contacts, quitting conbineonefitn");
                        goto quitrout;
                   }
                   badcont[ badcont[0] ] = j;
                }
                else if(laptype==2)
                   if(verbose>=1)
                      printf(" nterm of Fit %d  meets cterm of fit %d \n",ifit,j);
                else if(laptype==4)
                   if(verbose>=1)
                      printf(" Cterm of Fit %d meets nterm of fit %d \n",ifit,j);
#if 0
                else if(laptype==3)
                   printf("BAD FIT: nterm of Fit %d  meets nterm of fit %d \n",ifit,j);
                else if(laptype==5)
                   printf("BAD FIT:cterm of Fit %d  meets cterm  of fit %d\n",ifit,j);
                else
                   printf("FUNNY overlap:Fit %d res %d and fit %d res %d\n",ifit,ires,j,i2);
#endif
           }//nlap >0
           nextres:;
         }//j!=i
       }//j
   } // ires
   quitrout:;
}

void removebadfits(int iselect,struct griddata& den1,struct geometry& geo)
/* 
NOTE: THIS REMOVES ALL BAD FITS - MAY HAVE PROBLEM IF TOO MANY FITS OVERLAP
iselect = 0 - do all the checks
iselect = 1   just look for identical assigned fits
iselect = 2   look for identical assigned and unassigned fits
Looks for identical and bad fits and removes them
  TO DO - if fits do not overlap - should extend 1 so that it does
fuse[0][0,1] = data for nterm fusion
  fuse[0][0] = cterm end residue number of nterm fit that is going to be fused
  fuse[0][1] = nterm end residue number of cterm fit that is going to be fused
fuse[1][0,1] is same data for cterm fusion.
Fuse first at cterm end so that nterm data does not change
*/
{
   int i,ii,j,ncut;
   int ident[IDENTFITS][MAXIDENT][2],ibad[MAXIBAD][100][2],del[MAXDEL],iass;
   int iquit;

if(iselect ==1)
  goto identicalassign;
   printf("LOOKING FOR FITS THAT ARE IDENTICAL TO UNASSIGNED FITS:\n");
   del[0]=0; // keeps track of fits that should be deleted
   for(i=1;i<=fitnum[0];i++){ //First look for identical unassighed fits  - these must be short - not extended
     if(fit[i].val1>=0){ // not assigned
       iass= 0; //testing only unassigned versus all fits
       findidentical(&iquit,del,iass,i,ident);
       if(ident[0][0][0] == 1){ // 1 identical fit - expect them to be consecutive
         if(verbose>=1)
           printf(" Fit %d (residues %d to %d) overlaps fit %d in removebadfits1(residues %d to %d)\n",
             ident[1][0][0],ident[1][1][1],ident[1][ ident[1][0][1] ][1],
             i,ident[1][1][0],ident[1][ ident[1][0][1] ][0]);
         if( (ident[1][1][0] <= 1)&&
                 (ident[1][ ident[1][0][1] ][0] >=fit[i].nres-2) ){  // allow not ident at ends becaue end Ca can be in bad posiiton
            if(verbose>=1)
               printf("Fit %d overlaps all of fit %d - delete fit %d\n",ident[1][0][0],i,i);
              del[0]=del[0]+1;
              if(del[0]>= MAXDEL-1){
                 printf("ERROR - too many fits to delete - quitting combine\n");
                 goto quitrout;
              }
              del[ del[0] ] = i;
         }
         else{
           printf(" NOT COMPLETE OVERLAP - NEED TO EXTEND ONE FIT\n");
         }
       }
       else if (ident[0][0][0] >1){ // more than 1 identical fit - do not expect for unassigned fits
         printf("PROBLEM - %d identical fits to uassigned fit %d:\n",ident[0][0][0],i);
         for(j=1;j<=ident[0][0][0];j++)
            printf("  Fit %d (residues %d to %d)\n",
                ident[j][0][0],ident[j][1][1],ident[j][ ident[j][0][1] ][1]);
       }
     } // not assigned
     nexti:;
     if(iquit==1) // indicates array filled
          goto quitident;
   }// i
   quitident:;
   removenunordfits(fit,del);
   if(iselect ==2)
      goto identicalassign;
   printf("LOOKING FOR BAD UNASSIGNED FITS:\n");
   del[0]=0; // keeps track of fits that should be deleted
   for(i=1;i<=fitnum[0];i++){ //First look for bad unassighed fits  - these must be short - not extended
     if(fit[i].val1>=0){ // not assigned
       iass= 0; //testing only unassigned versus all fits
       findbad(&iquit,iass,i,ibad);
       if(ibad[0][0][0]>0){// if bad
           if(verbose>=1)
               printf(" Fit %d has bad contacts with %d fits:\n",i,ibad[0][0][0]);
           for(j=1;j<=ibad[0][0][0];j++){
              if(verbose>=1){              
                 printf("    Fit %d  has %d bad contacts with fit %d:  Residues ",
                          ibad[j][0][0],ibad[j][0][1],i);
                 for(ii=1;ii<=ibad[j][0][1];ii++)
                    printf(" (%d %d )",ibad[j][ii][1],ibad[j][ii][0]);
                 printf("\n");
              }
              if( (fit[ ibad[j][0][0] ].val1<0)||(ibad[0][0][0]>1) ){// if fit that it has bad contact with is assigned or has more than 1 bad fit, then delete
                 del[0]=del[0]+1;
                 if(del[0]>= MAXDEL-1){
                    printf("ERROR - too many fits to delete - quitting combine\n");
                    goto quitrout;
                 }
                 del[ del[0] ] = i;
                 if(verbose>=1)
                    printf(" DELETING fit %d because bad contacts\n",i);
                 goto nexti2;
              }
              else {//if reach here - have 1 bad contact with unassigned fit
                printf("WARNING: uassigned fits %d %d have bad contacts\n",i,ibad[1][0][0]);
              }
           }//j
        }//if bad
        if(iquit==1)
          goto quitbad;
     }// not assigned
     nexti2:;
  }// i
  quitbad:;
  removenunordfits(fit,del);
identicalassign:;
int ass1,ass2,end1,end2;
int fuse[2][2],ifuse[2][2],iend;
struct pepfit fit1;
float maxdist,mindist;
   fit1.nalloc = -10;
   allocatefit(fit1,30); // default allocation
   del[0]=0; // keeps track of fits that should be deleted
   printf("\n\n LOOKING FOR IDENTICAL TO ASSIGNED FITS:\n");
   iass=1; // indicates looking at assigned fits
   for(i=1;i<=fitnum[0];i++){ //After removing bad unassigned fits, look at assigned fits
//   for(i=4;i<=4;i++){ //After removing bad unassigned fits, look at assigned fits
     for(ii=1;ii<=del[0];ii++)
         if(del[ii]==i)
            goto nexti3;
     if(fit[i].val1<0){ //  assigned
       fuse[0][0] =fuse[1][0] = -1;  //indicates no overlap for fusion
       findidentical(&iquit,del,iass,i,ident);
       if(ident[0][0][0] >= 1){ // 1 identical fit - expect them to be consecutive
         if(verbose>=1)
             printf("Fit %d has identical regions with %d other assigned fits\n",i,ident[0][0][0]);
         for(ii=1;ii<=ident[0][0][0];ii++){
            if(verbose>=1)
               printf(" Fit %d (residues %d to %d) overlaps fit %d in removebadfits2(residues %d to %d)\n",
                ident[ii][0][0],ident[ii][1][1],ident[ii][ ident[ii][0][1] ][1],
                i,ident[ii][1][0],ident[ii][ ident[ii][0][1] ][0]);
            ass1 = -fit[ident[ii][0][0] ].val1 + ident[ii][1][1];
            ass2 = -fit[i].val1 + ident[ii][1][0];
            if(ass1 !=ass2){
               printf(" ERROR:  INCORRECT ASSIGNMENT FOR EITHER FIT %d or FIT %d\n",i,ident[ii][0][0]);
               printf("     THESE OVERLAPING FITS HAVE BEEN NOT BEEN FUSED\n");
               printf("     Sequence assignment for fit %d = %d; for fit %d = %d\n",ident[ii][0][0],ass1,i,ass2);
               printf("         These should be identical\n");
               fprintf(fperrmes," ERROR:  INCORRECT ASSIGNMENT FOR EITHER FIT %d or FIT %d\n",i,ident[ii][0][0]);
               fprintf(fperrmes,"     THESE OVERLAPING FITS HAVE BEEN NOT BEEN FUSED\n");
               fprintf(fperrmes,"     Sequence assignment for fit %d = %d; for fit %d = %d\n",ident[ii][0][0],ass1,i,ass2);
               fprintf(fperrmes,"         These should be identical\n");
               goto nexti3;
            }
            end1 = fit[i].nres-1 - ident[ii][ ident[ii][0][1] ][0];//The number of non-olverlpaed residue at end
            end2 = fit[ ident[ii][0][0] ].nres-1 -ident[ii][ ident[ii][0][1] ][1];
            if( (end1>=end2) && (ident[ii][1][0]>=ident[ii][1][1]) ){
                 if(verbose>=1)
                    printf("Fit %d overlaps all of fit %d - delete fit %d\n",
                        i,ident[ii][0][0],ident[ii][0][0]);
                 del[0]=del[0]+1;
                 if(del[0]>= MAXDEL-1){
                    printf("ERROR - too many fits to delete - quitting combine\n");
                    goto quitrout;
                 }
                 del[ del[0] ] = ident[ii][0][0];
                 goto nexti3; // this jumps out of loop - skips the possible fusion if overlap with another fit
            }
            else if( (end2>=end1)&& (ident[ii][1][1]>=ident[ii][1][0]) ){
                if(verbose>=1)
                   printf("Fit %d overlaps all of fit %d - delete fit %d\n",ident[ii][0][0],i,i);
                del[0]=del[0]+1;
                if(del[0]>= MAXDEL-1){
                   printf("ERROR - too many fits to delete - quitting combine\n");
                   goto quitrout;
                }
                del[ del[0] ] = i;
                goto nexti3; // this jumps out of loop - skips the possible fusion if overlap with another fit
           }
           else{// Not complete overlap - need to fuse - Ends of fit may be bad so fuse at residue 3 if possible
               if(verbose>=1)
                  printf(" NOT COMPLETE OVERLAP - FUSE THESE 2 FITS\n");
               if( ident[ii][1][0] <=3 ){ // then overlap at nterm end of i
                   ifuse[0][0]= ident[ii][0][0]; // nterm fit for fusion at nterm ofi
                   ifuse[0][1] = i;// cterm fit for fusion at nterm of i
                   fuse[0][1] = 1+(ident[ii][1][0]+ident[ii][ ident[ii][0][1] ][0])/2; // cterm fit fusion position
                   ncut = ident[ii][ ident[ii][0][1] ][0]-ident[ii][1][0] -fuse[0][1];
                   fuse[0][0] =  ident[ii][ ident[ii][0][1] ][1] -ncut; // nterm fit fusion position               }   
              }   
               if( ident[ii][ ident[ii][0][1] ][0] >= fit[i].nres-4 ){ // then overlap at cterm end of i
                   ifuse[1][0]= i; // nterm fit for fusion at cterm ofi
                   ifuse[1][1] = ident[ii][0][0];// cterm fit for fusion at cterm of i
                   fuse[1][1] = 1+(ident[ii][1][1]+ident[ii][ ident[ii][0][1] ][1])/2; // fusion position
                   ncut = ident[ii][ ident[ii][0][1] ][1]-ident[ii][1][1] -fuse[1][1]-1;
                   fuse[1][0] =  ident[ii][ ident[ii][0][1] ][0] - ncut;
                }   
            }
         } // ii
         if(fuse[1][0]>=0){ // then fuse at cterm end
             iend = 1; // ctermm end
             allocatefit(fit1,fit[ ifuse[iend][0]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][0] ],fit1,0,fuse[iend][0]);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][0] ],0,fit1.nres); // copy back to original fit
             allocatefit(fit1,fit[ ifuse[iend][1]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][1] ],fit1,fuse[iend][1]-1,fit[ ifuse[iend][1] ].nres);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][1] ],0,fit1.nres); // copy back to original fit
             geo.fitn = ifuse[iend][0];
             fusedist(fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ],&maxdist,&mindist);
             if(verbose>=2)
                printf("maxdist for fusion = %5.3f\n",maxdist);
             if(maxdist>geo.maxfusedist){
                printf("ERROR: CANNOT FUSE FITS %d and %d because fusedist = %5.3f is too large\n",
                        ifuse[iend][0],ifuse[iend][1],maxdist);
                goto nexti3;
             }
             fuseextfit2(den1,geo,0,fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ], ifuse[iend][1] );//"0" means remove 0 res from end of bestfit
             if(geo.quit==2){
                printf("ERROR: geo.quit=2 after fuseextfit2 - quitting routine\n");
                geo.quit=1;
                goto quitrout;
             }
             del[0]=del[0]+1;
             if(del[0]>= MAXDEL-1){
                    printf("ERROR - too many fits to delete - quitting combine\n");
                    goto quitrout;
             }
             del[ del[0] ] = ifuse[1][0];
             if(fuse[0][0]>=0){ // then will also fuse at nterm - must change fit number (fuseextfit() deletes lower numbered fit
                ifuse[0][1]= ifuse[1][1]; //the number that fuseextfit2 saved fused fit
             } 
         }
         if(fuse[0][0]>=0){ // then fuse at nterm end
             iend = 0; // nterm end
             allocatefit(fit1,fit[ ifuse[iend][0]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][0] ],fit1,0,fuse[iend][0]+1);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][0] ],0,fit1.nres); // copy back to original fit
             allocatefit(fit1,fit[ ifuse[iend][1]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][1] ],fit1,fuse[iend][1]-1,fit[ ifuse[iend][1] ].nres);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][1] ],0,fit1.nres); // copy back to original fit
             geo.fitn = ifuse[iend][0];
             fusedist(fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ],&maxdist,&mindist);
             if(verbose>=2)
                printf("maxdist for fusion = %5.3f\n",maxdist);
             if(maxdist>geo.maxfusedist){
                printf("ERROR: CANNOT FUSE FITS %d and %d because fusedist = %5.3f is too large\n",
                        ifuse[iend][0],ifuse[iend][1],maxdist);
                goto nexti3;
             }
             fuseextfit2(den1,geo,0,fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ], ifuse[iend][1] );//"0" means remove 0 res from end of bestfit
             if(geo.quit==2){
                printf("ERROR: geo.quit=2 after fuseextfit2 - quitting routine\n");
                geo.quit=1;
                 goto quitrout;
             }
             del[0]=del[0]+1;
             if(del[0]>= MAXDEL-1){
                    printf("ERROR - too many fits to delete - quitting combine\n");
                    del[0]=del[0]-1;
                    goto quitrout;
             }
             del[ del[0] ] = ifuse[0][1];
         }
       }// >=1
     }//if assigned
     nexti3:;
     if(iquit==1) 
        goto quitrout;
   }//i
  quitrout:;
  removenunordfits(fit,del);
  freepepfit(fit1);
}//removebadfits()


void removeonebadfit(int ifit,int iselect,struct griddata& den1,
       struct geometry& geo,int del[MAXDEL],int origfit)
/* 
origfit = fit # that new ifit was made from bf by symop
Similar to removebadfit execpt just checks ifit, not all the fits - and
  ASSUMES BOTH FITS ARE ASSIGNED.
Fuse first at cterm end so that nterm data does not change
*/
{
   int i,ii,id;
   int ident[IDENTFITS][MAXIDENT][2],iass;
   int iquit,initdel;
int ass1,ass2,end1,end2;
int fuse[2][2],ifuse[2][2],iend;
struct pepfit fit1;
float maxdist,mindist;


 i = ifit;
   fit1.nalloc = -10;
   initdel = del[0];
   allocatefit(fit1,30); // default allocation
   iass=0; // NEW - test all fits - assigned and unassigned fits
//   iass=1; // indicates looking at assigned fits
//     if(fit[i].val1<0){ // NEW - test all fits - assigned and unassighed 
       fuse[0][0] =fuse[1][0] = -1;  //indicates no overlap for fusion
//if(i==130)
//  savefit();
       findoneidentical2(&iquit,del,iass,i,ident);
       if(ident[0][0][0] >= 1){ // 1 or more  identical fit - expect them to be consecutive
         if(verbose>=1){
           printf("Fit %d has identical regions with %d other assigned fits:",i,ident[0][0][0]);
           for(ii=1;ii<=ident[0][0][0];ii++)
              printf(" %d",ident[ii][0][0]);
           printf("\n");
         }
         if(fit[i].val1>=0){ // if fit i is not assighed
             for(ii=1;ii<=ident[0][0][0];ii++){ //first check if overlaps assigned fit
               for(id=1;id<=del[0];id++) // marked for deletion
                  if(del[id]== ident[ii][0][0])
                    goto nextii1;
               if(fit[ ident[ii][0][0] ].val1<0){
                   if(verbose>=1)
                      printf(" Assigned Fit %d (residues %d to %d) overlaps unassigned fit %d (residues %d to %d)\n",
                       ident[ii][0][0],ident[ii][1][1],ident[ii][ ident[ii][0][1] ][1],
                       i,ident[ii][1][0],ident[ii][ ident[ii][0][1] ][0]);
                   if(verbose>=1)
                      printf("Delete unassigned fit %d\n",i);
                   del[0]=initdel+1; // only delete fit i
                   del[ del[0] ] = i;
                   goto quitrout; // this jumps out of loop - i is completely overlaped and deleted
               }
               nextii1:;
             }//ii
             for(ii=1;ii<=ident[0][0][0];ii++){ //then check if overlaps unassigned fit
               for(id=1;id<=del[0];id++) // marked for deletion
                  if(del[id]== ident[ii][0][0])
                    goto nextii2;
               if(fit[ ident[ii][0][0] ].val1>=0){ // At this point this should always be true
                   if(verbose>=1)
                      printf(" Unassigned  unassigned Fit %d (residues %d to %d) overlaps unassigned fit %d (residues %d to %d)\n",
                          ident[ii][0][0],ident[ii][1][1],ident[ii][ ident[ii][0][1] ][1],
                          i,ident[ii][1][0],ident[ii][ ident[ii][0][1] ][0]);
                   if(fit[ ident[ii][0][0] ].nres > fit[i].nres){
                        if(verbose>=1)
                           printf("Delete  unassigned fit %d because shorter than fit %d\n",i,ident[ii][0][0]);
                        del[0]=del[0]+1; // only delete fit i
                        del[ del[0] ] = i;
                        goto quitrout; // this jumps out of loop - i is completely overlaped and deleted
                   }
                   else{
                        if(verbose>=1)
                            printf("Delete fit %d because shorter than fit %d\n",ident[ii][0][0],i);
                        del[0]=del[0]+1; // only delete fit i
                        del[ del[0] ] = ident[ii][0][0];
                    }
               }// val1>=0
               nextii2:;
             }//ii
             goto quitrout;  // If i unassgined, quit
         } // i not assigned
         for(ii=1;ii<=ident[0][0][0];ii++){ // Use this section for assigned fit i
            for(id=1;id<=del[0];id++) // marked for deletion
               if(del[id]== ident[ii][0][0])
                  goto nextii;
            if(verbose>=1)
              printf(" Fit %d (residues %d to %d) overlaps fit %d in removeonebadfit (residues %d to %d)\n",
                ident[ii][0][0],ident[ii][1][1],ident[ii][ ident[ii][0][1] ][1],
                i,ident[ii][1][0],ident[ii][ ident[ii][0][1] ][0]);
//if( (ident[ii][0][0]==77)&&(i==130) )
// printf("\n");
            if(fit[ident[ii][0][0] ].val1>=0){
                 if(verbose>=1)
                     printf("Delete unassigned fit %d because overlapped by assigned fit %d\n",ident[ii][0][0],i); 
                 del[0]= del[0]+1;
                 del[ del[0] ] = ident[ii][0][0];
                 goto nextii;
            }
            ass1 = -fit[ident[ii][0][0] ].val1 + ident[ii][1][1];
            ass2 = -fit[i].val1 + ident[ii][1][0];
            if(ass1 !=ass2){
               printf(" ERROR:  INCORRECT ASSIGNMENT FOR EITHER FIT %d or FIT %d - FIT %d has been deleted\n",i,ident[ii][0][0],i);
               printf("     Sequence assignment for fit %d is %d; For fit %d (made by rot/tran on fit %d)= %d\n",ident[ii][0][0],ass1,i,origfit,ass2);
               printf("         These should be identical\n");
               fprintf(fperrmes," ERROR:  INCORRECT ASSIGNMENT FOR EITHER ROTATED FIT %d or FIT %d \n",origfit,ident[ii][0][0]);
               fprintf(fperrmes,"     Sequence assignment for fit %d is %d; For fit %d = %d\n",ident[ii][0][0],ass1,origfit,ass2);
               fprintf(fperrmes,"         These should be identical\n");
               del[0]=initdel+1; // only delete fit i- do not delete any other fit
               del[ del[0] ] = i;
               goto quitrout; // this jumps out of loop - i is completely overlaped and deleted
            }
            end1 = fit[i].nres-1 - ident[ii][ ident[ii][0][1] ][0];//The number of non-olverlpaed residue at end
            end2 = fit[ ident[ii][0][0] ].nres-1 -ident[ii][ ident[ii][0][1] ][1];
            if( (end1>=end2) && (ident[ii][1][0]>=ident[ii][1][1]) ){
                 if(verbose>=1)
                    printf("Fit %d overlaps all of fit %d - delete fit %d\n",
                        i,ident[ii][0][0],ident[ii][0][0]);
                 del[0]= del[0]+1;
                 del[ del[0] ] = ident[ii][0][0];
            }
            else if( (end2>=end1)&& (ident[ii][1][1]>=ident[ii][1][0]) ){ // fit i is completely overlapped
                if(verbose>=1)
                   printf("Fit %d overlaps all of fit %d - delete fit %d\n",ident[ii][0][0],i,i);
                del[0]=initdel+1; // only delete fit i- do not delete any other fit
                del[ del[0] ] = i;
                goto quitrout; // this jumps out of loop - i is completely overlaped and deleted
           }
           else{// Not complete overlap - need to fuse - Ends of fit may be bad so fuse at residue 3 if possible
               if( (fit[i].nres-ident[ii][1][0])>(fit[ ident[ii][0][0] ].nres-ident[ii][1][1]) ){ // then overlap at nterm end of i
//              if( ident[ii][1][0] <=3 ){ // then overlap at nterm end of i
                   if(verbose>=1)
                       printf(" NOT COMPLETE OVERLAP - FUSE THESE 2 FITS - fit %d is cterm fit\n",i);
                   ifuse[0][0]= ident[ii][0][0]; // nterm fit for fusion at nterm ofi
                   ifuse[0][1] = i;// cterm fit for fusion at nterm of i
                   if( ident[ii][0][1]==1)// if only 1 resdidue overlap - then fuse at it
                       fuse[0][1]=ident[ii][ ident[ii][0][1] ][0];
                   else
                       fuse[0][1] =(ident[ii][1][0]+ident[ii][ ident[ii][0][1] ][0])/2; // cterm fit fusion position
                   fuse[0][0] =ident[ii][1][1] + (fuse[0][1]-ident[ii][1][0]); // residue in fit i corresponding to fuse[1][1]
//                   ncut = ident[ii][ ident[ii][0][1] ][0]-ident[ii][1][0] -fuse[0][1];
//                   fuse[0][0] =  ident[ii][ ident[ii][0][1] ][1] -ncut; // nterm fit fusion position               }   
               }
               else{
//               if( ident[ii][ ident[ii][0][1] ][0] >= fit[i].nres-4 ){ // then overlap at cterm end of i
                   if(verbose>=1)
                       printf(" NOT COMPLETE OVERLAP - FUSE THESE 2 FITS - fit %d is nterm fit\n",i);
                   ifuse[1][0]= i; // nterm fit for fusion at cterm ofi
                   ifuse[1][1] = ident[ii][0][0];// cterm fit for fusion at cterm of i
                   if( ident[ii][0][1]==1)// if only 1 resdidue overlap - then fuse at it
                       fuse[1][1]=ident[ii][ ident[ii][0][1] ][1];
                   else // fuse at intermediate residue
                      fuse[1][1] = 1+(ident[ii][1][1]+ident[ii][ ident[ii][0][1] ][1])/2; //overlap fuse[1][1]-1 to fuse[1][1] of cterm fitfusion position
                   fuse[1][0] =1+ ident[ii][1][0] + (fuse[1][1]-ident[ii][1][1]); // residue in fit i corresponding to fuse[1][1]
//                   ncut = ident[ii][ ident[ii][0][1] ][1]-ident[ii][1][1] -fuse[1][1]-1;
//                   fuse[1][0] =  ident[ii][ ident[ii][0][1] ][0] - ncut;
               }   
            }
            nextii:;
         } // ii
         if(fuse[1][0]>=0){ // then fuse at cterm end
             iend = 1; // ctermm end
             allocatefit(fit1,fit[ ifuse[iend][0]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][0] ],fit1,0,fuse[iend][0]);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][0] ],0,fit1.nres); // copy back to original fit
             allocatefit(fit1,fit[ ifuse[iend][1]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][1] ],fit1,fuse[iend][1]-1,fit[ ifuse[iend][1] ].nres);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][1] ],0,fit1.nres); // copy back to original fit
             geo.fitn = ifuse[iend][0];
             fusedist(fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ],&maxdist,&mindist);
             if(verbose>=2)
                printf("maxdist for fusion = %5.3f\n",maxdist);
             if(maxdist>geo.maxfusedist){
                printf("ERROR: CANNOT FUSE FITS %d and %d because fusedist = %5.3f is too large\n",
                        ifuse[iend][0],ifuse[iend][1],maxdist);
                goto quitrout;
             }
             fuseextfit2(den1,geo,0,fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ], ifuse[iend][1] );//"0" means remove 0 res from end of bestfit
             del[0]=del[0]+1;
             if(del[0]>= MAXDEL-1){
                    printf("ERROR - too many fits to delete - quitting combine\n");
                    goto quitrout;
             }
             if(ifuse[1][0]<ifuse[1][1]){
                del[ del[0] ] = ifuse[1][0];
             }
             else{ 
                del[ del[0] ] = ifuse[1][1];
             } 
   savefit(); // helps in debugging
         }
         if(fuse[0][0]>=0){ // then fuse at nterm end - Even if already fused at cterm - this should work
             iend = 0; // nterm end
             allocatefit(fit1,fit[ ifuse[iend][0]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][0] ],fit1,0,fuse[iend][0]+1);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][0] ],0,fit1.nres); // copy back to original fit
             allocatefit(fit1,fit[ ifuse[iend][1]  ].nres +2); // make sure big enough allocation
             copyfittofit(fit[ ifuse[iend][1] ],fit1,fuse[iend][1]-1,fit[ ifuse[iend][1] ].nres);// make truncated nterm end
             copyfittofit(fit1,fit[ ifuse[iend][1] ],0,fit1.nres); // copy back to original fit
             geo.fitn = ifuse[iend][0];
             fusedist(fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ],&maxdist,&mindist);
             if(verbose>=2)
                printf("maxdist for fusion = %5.3f\n",maxdist);
             if(maxdist>geo.maxfusedist){
                printf("ERROR: CANNOT FUSE FITS %d and %d because fusedist = %5.3f is too large\n",
                        ifuse[iend][0],ifuse[iend][1],maxdist);
                goto quitrout;
             }
             fuseextfit2(den1,geo,0,fit[ ifuse[iend][0] ],fit[ ifuse[iend][1] ], ifuse[iend][1] );//"0" means remove 0 res from end of bestfit
             del[0]=del[0]+1;
             if(del[0]>= MAXDEL-1){
                    printf("ERROR - too many fits to delete - quitting combine\n");
                    goto quitrout;
             }
             if(ifuse[0][1]<ifuse[0][0])
                del[ del[0] ] = ifuse[0][1];
             else 
                del[ del[0] ] = ifuse[0][0];
   savefit(); // helps in debugging
         }
       }// >=1
//     }//if assigned
  quitrout:;
  freepepfit(fit1);
}


void combine(struct griddata& den1,struct geometry& geo,char combinefile[])
/* General program for combining severalfit files
*/
{
   int iselect;

  readcombinefile(geo,combinefile);
  iselect = 0; // indicates use all checks in removebad
  removebadfits(iselect,den1,geo);
  quitrout:;
}

void saveonesym(FILE *fp,float tran[],float rot[3][3])
// Prints to file the translation and rotation matrix
// Modified from printsymoperator
{
   fprintf(fp,"%5.3f %5.3f %5.3f\n",tran[0],tran[1],tran[2]);
   fprintf(fp,"%5.3f %5.3f %5.3f\n",rot[0][0],rot[0][1],rot[0][2]);    
   fprintf(fp,"%5.3f %5.3f %5.3f\n",rot[1][0],rot[1][1],rot[1][2]);    
   fprintf(fp,"%5.3f %5.3f %5.3f\n",rot[2][0],rot[2][1],rot[2][2]);    
}

void readonesym(FILE *fp,float tran[],float rot[3][3])
{
   char line[300];
   
   fgets(line,200,fp); /* read each line*/
   sscanf(line,"%f %f %f",&tran[0],&tran[1],&tran[2]);
   fgets(line,200,fp); /* read each line*/
   sscanf(line,"%f %f %f\n",&rot[0][0],&rot[0][1],&rot[0][2]);    
   fgets(line,200,fp); /* read each line*/
   sscanf(line,"%f %f %f\n",&rot[1][0],&rot[1][1],&rot[1][2]);    
   fgets(line,200,fp); /* read each line*/
   sscanf(line,"%f %f %f\n",&rot[2][0],&rot[2][1],&rot[2][2]);    
}

void savensym(char symfile[],int nsym,float rot[NSYM][3][3],float tran[NSYM][3] )
{
   int i;
   FILE *fp;

   fp= fopen(symfile,"w");
   fprintf(fp,"%d\n",nsym);
   for(i=1;i<=nsym;i++){
     fprintf(fp,"Rotation/Translation #%d:\n",i);
     saveonesym(fp,tran[i],rot[i]);
   }
   fclose(fp);
}
   
void readnsym(char symfile[],int *nsym,float rot[NSYM][3][3],float tran[NSYM][3] )
{
   int i;
   char line[300];
   FILE *fp;

   printf("READING SYMFILE  %s\n",symfile);
   if( (fp=fopen(symfile,"r") ) ==NULL) { // check if file exists
         printf("ERROR: The %s file does not exist",symfile);
   }
   else{
     fgets(line,200,fp); /* read each line*/
     sscanf(line,"%d",nsym);
     printf("  Number of sym operators = %d\n",*nsym);
     for(i=1;i<= *nsym;i++){
        fgets(line,200,fp); //read and discard verbiage
        readonesym(fp,tran[i],rot[i]);
     }
   }
   fclose(fp);
}


void expand(struct griddata& den1,struct geometry& geo, const char expandfile[])
/*
NEW input fitfile= expandfile instead of combinefile
 General program expanding a fit file into sym. related regions.
*/
{
   int i,isym[NSYM][2],icheck,ifind,iselect,iplus;
   int nsym,del[MAXDEL];
   float brot[NSYM][3][3],btran[NSYM][3];
   char savfile[300];
   char errmesfile[300];
   float maxfusedist;

  maxfusedist = geo.maxfusedist;
  geo.maxfusedist = 1.8;
  printf("STARTING EXPAND USING geo.maxfusedist = %5.3f\n",geo.maxfusedist);
  strcpy(errmesfile,maidfile.name);
  strcat(errmesfile,".expanderror");
  if( (fperrmes=fopen(errmesfile,"w") ) ==NULL) {
      fprintf(stderr,"error opening file %s\n",errmesfile);
      exit(0);
  }
  icheck = 0; // Check all unassigned fits (=ilim = 25 residues) of assigned fits
//  readcombinefile(geo,expandfile);
  readfit(geo,expandfile);
  iselect = 0; // indicates use all checks in removebad
  nsym=0; // Number of rot/tran operators found
  iplus = 0; // Means try all ufits - do not check if rot/tran unique
//  iplus = 1; // Means try only  ufits >i, prevents duplication

  removebadfits(iselect,den1,geo); // Look at all initial fits for bad or overlaping fits before expanding
#if 0  // There should not be any sym. between assigned and unassgined, since assigned should have been used to assign the unassinged
  printf("\n\n LOOKING FOR SYMMETRY BETWEEM ASSIGNED AND UNASSIGNED FITS:\n");
//  iplus = 1; // Means try only  ufits >i, prevents duplication
  iplus = 0; // Means try all ufits - do not check if rot/tran unique
  for(i=1;i<=fitnum[0];i++){ //Look for symetry between assigned fits and unassigned fits 
     if(fit[i].val1<0){ // assigned
         ifind =1; // return best rot matrix after finding first sym mate
         symfit(iplus,&ifind,icheck,geo,den1,i,isym,&nsym,brot,btran,del);
         if(nsym==NSYM-1){
            printf("ERROR: After symfit -Reached allowed limited of number of rot/tran operators\n");
            printf("    Increase NSYM and recompile\n");
           goto quitrout;
         }
     }
  }
#endif
  printf("\n\n LOOKING FOR SYMMETRY BETWEEM ASSIGNED  FITS:\n");
  del[0]=0; // keeps track of fits that should be deleted
  for(i=1;i<=fitnum[0];i++){ //Look for symetry between assigned fits 
     if(fit[i].val1<0){ // assigned
         ifind =1; // return best rot matrix after finding first sym mate
         symfitass(iplus,&ifind,icheck,geo,den1,i,isym,&nsym,brot,btran,del);
         if(nsym==NSYM-1){
            printf("ERROR: After symfitass -Reached allowed limited of number of rot/tran operators\n");
            printf("    Increase NSYM and recompile\n");
           goto quitrout;
         }
     }
  }
  if(verbose>=1){
     printf(" Saved %d rot/tran fits in file %s\n",nsym,maidfile.nsym);
     savensym(maidfile.nsym,nsym,brot,btran);
  }  
//  del[0]=0; // eset to 0 keeps track of fits that should be deleted
  quitrout:;
  applyrottran(geo,den1,nsym,brot,btran,del); // Apply rot/tran opertors to each fit, save if in good density - fuse and find del[]
  removenunordfits(fit,del);
  printf("FINISHED EXPAND\n");
  strcpy(savfile,maidfile.fitfile);
  strcat(savfile,"_expand");
  savebackupfit(savfile); // copies fitfile to savfile      
  printf("LOOK IN FILE %s FOR SUMMARY OF PROGLEM ASSIGNMENTS\n",errmesfile);
  geo.maxfusedist=maxfusedist;
  fclose(fperrmes);
}

void findfrsphere(struct pepfit fit1,float *frsphere)
/* Adapted from saveonefit
*/
{
   int j,ii,jj,is,count,totcount,rottype,acount,atomname;
   float dist;

    count = totcount = 0;
    for(j=0;j<fit1.nres;j++){ // each residue
       if( (j==0)||(j==fit1.nres-1)||(fit1.nside[j]>=0) ){// then side not residue
         for(ii=1;ii<=4+fit1.nside[j];ii++){ // all atom positions per reisude 
           totcount = totcount+1;
           for(is =1;is<=maidfile.nsphere;is++){
              dist = distv(fit1.pos[j][ii],sp[is].cenpos);
                if(dist <= sp[is].radius){
                   count = count+1;
                   goto nextii;
                }
           }// is
           nextii:;
         }//ii
       } // side not residue
       else{ // side = residue
          for(ii=1;ii<=4;ii++){ // main chain atoms
            totcount = totcount+1;
            for(is =1;is<=maidfile.nsphere;is++){
              dist = distv(fit1.pos[j][ii],sp[is].cenpos);
                if(dist <= sp[is].radius){
                   count = count+1;
                   goto nextii2;
                }
            }// is
            nextii2:;
          }//ii
          rottype=  -fit1.nside[j]; // residue name (number)
          acount = 0;
          for(ii=1;ii<=side[rottype][0][0]+1;ii++){ // copy each side chain branch
            for(jj=0;jj<side[rottype][ii][0];jj++){
              atomname=side[rottype][ii][jj+1];
              if(atomname>5){ // not part of mainchain or Cb
                acount = acount+1;
                totcount = totcount+1;
                for(is =1;is<=maidfile.nsphere;is++){
                  dist = distv(fit1.pos[j][acount+5],sp[is].cenpos);
                  if(dist <= sp[is].radius){
                     count = count+1;
                     goto nextjj;
                  }
                } // is
              } // atopmname
              else{ // = 5 = Cb
                totcount = totcount+1;
                for(is =1;is<=maidfile.nsphere;is++){
                  dist = distv(fit1.pos[j][atomname],sp[is].cenpos);
                  if(dist <= sp[is].radius){
                     count = count+1;
                     goto nextjj;
                  }
                } // is
              }// Cb
              nextjj:;
            }// jj
          }//ii
       } // side = residue
    } // j  each residue 
    *frsphere = (float)count/totcount;
}    

#if 0 // OLD version
void assignseqment(struct geometry& geo)
/*   assign fits to proper segment as follow:
    Define segment number = -fit[i].val2 + 1000  (i.e. val2 = -1001 = segment 1(
    For each fit:
       Store val1 and find frsphere - the fraction of atoms that lie in sphere.
         For each val1: 
            Assign seqment # 0 (= -val2) to the fit with larges frsphere.
         Then for other fits with same val1 -
            Apply all sym operators i = 1,2,.. and find frsphere -
              The assign the segment # = i value with larges frsphere.
  sphere data: i = 1,... maidfile.nsphere
     sp[i].radius; sp[i].cenpos[0,1,2]
*/
{
   int i,j,imax,ifit[NSYM+5],count,nsym,isym,maxside;
   float frsphere,maxsphere,rot[NSYM][3][3],tran[NSYM][3];
   struct pepfit fit2;

   fit2.nalloc = -10;
   allocatefit(fit2,fit[1].nres +2); // make sure big enough allocation
   maxside = 3; // max. number of side chain atoms to use if unassighed.
   recallspheresbat(geo);
   printf(" STARTING ROUTINE assignseqment()\n");
   readnsym(maidfile.nsym,&nsym,rot,tran); // read rot/tran operators
   for(i=1;i<=fitnum[0];i++){
     if( (fit[i].val1 < 0)&&(fit[i].val2 > -999) ){ // sequence assigned and segment not assigned
        count = 1;
        findfrsphere(fit[i],&frsphere);
        if(verbose>=2)
            printf("frpshere[%d] = %5.3f\n",i,frsphere);
        ifit[count]=i;
        maxsphere = frsphere;
        imax = ifit[count];
        for(j=i+1;j<=fitnum[0];j++){
          if( (fit[j].val1 < 0)&&(fit[j].val2 > -999) ){
            if( ( (-fit[j].val1 <= -fit[i].val1)&&
                 (-fit[j].val1+fit[j].nres-1>= -fit[i].val1) ) ||
            ( (-fit[i].val1 <= -fit[j].val1)&&
                 (-fit[i].val1+fit[i].nres-1>= -fit[j].val1) ) ) {//  then fit i and j have similar seq. assignment
              count=count+1;
              ifit[count] = j;               
              findfrsphere(fit[j],&frsphere);
              if(verbose>=2)
                 printf("fit %d frpshere = %5.3f\n",j,frsphere);
              if(frsphere>maxsphere){ //Then make fit j the maximum 1
                  maxsphere = frsphere;
                  imax = ifit[count]; // = j
              }
            }
          }// if
        }//j
        fit[imax].val2 = -1000;  // segment = 0
        if(nsym< count-1){
            printf("ERROR: THERE ARE FEWER SYM. OPERATORS THAN FITS WITH SAME SEQUNECE ASSIGNMENT\n");
        }
        for(isym = 1;isym <=nsym;isym++){  // find operator that give largest frsym 
           maxsphere = -100.0; 
           for(j=1;j<=count;j++){ // for rest - assign segment # = sym oper. that produces largest frsphere
              if( (ifit[j]!=imax)&&(fit[ ifit[j] ].val2 > -999) ){
                 movefitaa(fit[ ifit[j] ],fit2,rot[isym],tran[isym],maxside);
                 findfrsphere(fit2,&frsphere);
                 if(verbose>=1)
                     printf("rotated fit %d frpshere = %5.3f\n",ifit[j],frsphere);
                 if(frsphere>maxsphere)
                     imax = ifit[j];
              }//ifit
            }// j
            fit[imax].val2 = -1000 - isym;
        } // isym
      }// seq. assigned
      else if(fit[i].val2 > -999){  // sequence and segment not assigned - assignme val2  = -2000
        fit[i].val2 = -2000;
     } // unassigned 
  }
  savefit(); //save the val2 segment assignments          
  freepepfit(fit2);
  printf("FINISHED ROUTINE assignseqment()\n");
}
#endif


#if 1 // NEW version
void assignseqment(struct griddata& den1,struct geometry& geo, const char callfile[])
/* 
25jan01:  callfile is the process that called this - either extend or expand
    This name is included in pdb file name
5jan01:  First look for bone connections to other fits.
16nov00:  STARTED NEW PROCEDURE FOR ASSIGNING SEGMENT:
   1) Find longest fit.  (If more than one fit with this length choose the one
            with the largest fraction in sphere.
        Call this tfit = 0, and assign it seg = 0;
   2) Go through all assigned fits and find the fit that: 1)does not overlap seq=0 and
         2) end Ca is closest to either end of seq=0 - and assign it seg=0
   3) Continue process in 2 for all fits (using all the ends in seq=0)
   4) Apply all the sym operators to this segment =0 
         assign the seg = count to any fit (assigned and unassigned) that it overlaps  where 
         count NCS operators that overlap some fit.
   5) Hopefully, this will have assigned segment to all fits that have sequence assigned,
        Check - if finds assigned sequence with segment not assigned - assign it seg.,
        and apply NCS opeators - and continue until cannot find anymore.
   6)Then check all fits with sequence and segment not assinged (this means cannot overlap
          previously assigned fit) and repeat 2.

  sphere data: i = 1,... maidfile.nsphere
     sp[i].radius; sp[i].cenpos[0,1,2]
*/
{
   int i,j,jj,ii,maxn,imax[50],maxside,amax,nsym,iass;
   int count,del[MAXDEL],iquit;
   float frsphere,maxfr,brot[3][3],btran[3];;
   int ident[IDENTFITS][MAXIDENT][2],addone,origfitnum;
   float rotm[3][3],tran[3],frout,finalden,minden,minden0,mindist,dist;
   int imin,nterm,nterm2,cterm,cterm2,nseg[NEND]; // nseg[0] =number of fits assigned to count = 1, nseg[i] = fit number
   int maxres,minres,foundcon,icon;
   char fitpdbfile[300];
   int storecon[20]; //New 26aug00 - prevents endless loop
   FILE *fp;

   del[0]=0; // Note this is not used in this routine, but needed by subroutine that are called keeps track of fits that should be deleted
   minden0 = 20.0;  // Minimum average density required for potential sym. fit
   maxside = 3; // max. number of side chain atoms to use if unassighed.
   readfit(geo,maidfile.fitfile);
   origfitnum = fitnum[0];
   recallspheresbat(geo);
   printf(" STARTING ROUTINE assignseqment()\n");
//   readnsym(maidfile.nsym,&nsym,brot,btran); // read rot/tran operators
   maxn = -100;
   imax[0]=0;
   for(i=1;i<=origfitnum;i++){ // find longest fit with uassaigned segment
       if(fit[i].val2 > -900){ // segment not yet assigned
          if(fit[i].nres>maxn){
              maxn = fit[i].nres;
              imax[0]= 1; // counts number of identical with same max length
              imax[1] = i;
          }
          else if (fit[i].nres==maxn){
             imax[0]=imax[0]+1; 
             imax[ imax[0] ] = i;
             if(imax[0]>48){
               printf("ERROR: TOO many equal lenght fits in assignseqment()\n");
               goto quitrout;
            } 
          }
        }// if > -999
   }// i
   if(imax[0]==0){ // all fits have segment assigned
        printf("All fits have segemnts already assigned\n");
        goto quitrout;
   }
   else if(imax[0]>1){ // more than 1 max - find the 1 with largest frsphre = amax
        maxfr = -10.0;
        for(j=1;j<=imax[0];j++){
           findfrsphere(fit[ imax[j] ],&frsphere);
           if(frsphere>maxfr){
              amax = imax[j];
              maxfr = frsphere;
           }
        }
   }
   else{ //  imax[0]=1
       amax = imax[1];
   }
   count = 1; // = segment number
   fit[amax].val2 = -900 - count;  // segment number = -(val2-900)
   nseg[0]=1; // number of fits in count = 0
   nseg[ nseg[0] ] = amax;
   geo.fixend=2;  // extend at cterm end
   foundcon =1;// indicates found connection
   icon = amax; // fit used to look for cterm or nterm connection
   storecon[0]=1;
   storecon[1]=icon;//Stores tested fit so do not get into endless loop
   while(foundcon>0){ // Look for connection to other fits leading from cterm end
      batchextendcomb(den1,geo,icon,&foundcon);
      if(foundcon>0){// assign connected fit to
          for (i=1;i<=storecon[0];i++)
            if (foundcon==storecon[i])
               goto quitwhile1; 
          storecon[0]=storecon[0]+1;
          if (storecon[0]>18){
              printf("ERROR: storecon[0]>18 - quitting assignsegment\n");
              goto quitwhile2;
          }
          storecon[ storecon[0] ]=foundcon;
          fit[foundcon].val2 = -900 - count;  // add this segment to seg 1
          nseg[0]=nseg[0]+1;
          nseg[ nseg[0] ] = foundcon;
          icon = foundcon; // continue try to extend using new fit 
      }   
   }
   quitwhile1:;
   geo.fixend=1;  //Repeat the process  extending at nterm end
   foundcon =1;// indicates found connection
   icon = amax; // fit used to look for cterm or nterm connection
   while(foundcon>0){ // Look for connection to other fits leading from cterm end
      batchextendcomb(den1,geo,icon,&foundcon);
      if(foundcon>0){// assign connected fit to 
          for (i=1;i<=storecon[0];i++)
            if (foundcon==storecon[i])
               goto quitwhile2; 
          storecon[0]=storecon[0]+1;
          if (storecon[0]>18){
              printf("ERROR: storecon[0]>18 - quitting assignsegment\n");
              goto quitwhile2;
          }
          storecon[ storecon[0] ]=foundcon;
          fit[foundcon].val2 = -900 - count;  // add this segment to seg 1
          nseg[0]=nseg[0]+1;
          nseg[ nseg[0] ] = foundcon;
          icon = foundcon; // continue try to extend using new fit 
      }   
   }
   quitwhile2:;
    while(1){ //Assign remaining non-overlaping assinged fits to count = 1 ( use goto to leave loop)
     mindist =  100000.0;
     for(i=1;i<=origfitnum;i++){ // find closest fit with uassaigned segment
        if( (fit[i].val2> -900)&&(fit[i].val1<0) ){ //sequence assigned and segment not assigned
          for(ii=1;ii<= nseg[0];ii++){ // check if overlaps seg=1
             nterm = -fit[ nseg[ii] ].val1;
             cterm = nterm + fit[ nseg[ii] ].nres-2;
             nterm2 = -fit[i].val1;
             cterm2 = nterm2+fit[i].nres-2;
             if( ((nterm2>=nterm)&&(nterm2<=cterm)) ||
                    ((cterm2>=nterm)&&(cterm2<=cterm)) ){// then overlaps
                 goto nexti;
             }
          }//ii
          for(ii=1;ii<= nseg[0];ii++){ // find fit with shortest gap distg distance
             dist= distv( fit[ nseg[ii] ].pos[0][1],fit[i].pos[ fit[i].nres-1 ][1] ); // cterm2,nterm dist
//printf(" at 1 i=%d ii=%d dist = %5.3f\n",i,ii,dist);
             if(dist<mindist){
               mindist = dist;
               imin = i;
             }
             dist= distv( fit[ nseg[ii] ].pos[ fit[nseg[ii]].nres-1 ][1],fit[i].pos[0][1] ); // nterm2,cterm dist
//printf(" at 2 i=%d  dist = %5.3f\n",i,dist);
             if(dist<mindist){
               mindist = dist;
               imin = i;
             }
           } //ii
        } // not assigned   
         nexti:;
      }//i  
      if(mindist < 1000.0){ // indicate found at another segment to add
              fit[imin].val2 = -900 - count;  // add this segment to seg 1
              nseg[0]=nseg[0]+1;
              nseg[ nseg[0] ] = imin;
       }
       else{ // no more fits that can be assigned
             goto quitwhile;
       }
   }// while         
   quitwhile:;
#if 0
   while(1){ //Assign closest unassigned fit to seg=1 ( use goto to leave loop)
     mindist =  100000.0;
     for(i=1;i<=origfitnum;i++){ // find closest fit with uassaigned segment
        if( (fit[i].val2> -900)&&(fit[i].val1>0)&&(fit[i].val1!=100000) ){ // segment not assigned (sequence also must be not assigned
          for(ii=1;ii<= nseg[0];ii++){ // find fit with shortest gap distg distance
             dist= distv( fit[ nseg[ii] ].pos[0][1],fit[i].pos[ fit[i].nres-1 ][1] ); // cterm2,nterm dist
             if(dist<mindist){
               mindist = dist;
               imin = i;
             }
             dist= distv( fit[ nseg[ii] ].pos[ fit[nseg[ii]].nres-1 ][1],fit[i].pos[0][1] ); // nterm2,cterm dist
             if((dist<mindist)&&(dist<20.0) ){ // require that unassigned fit must be closer than 20 Ang.
               mindist = dist;
               imin = i;
             }
           } //ii
           if(mindist < 20.0){ // indicate found  another segment to add
              fit[imin].val2 = -900 - count;  // add this segment to seg 1
              nseg[0]=nseg[0]+1;
              nseg[ nseg[0] ] = imin;
              refinerottran(geo,den1,brot[isym],btran[isym],maxside,rotm,tran,imin,
                  minden,&finalden,&frout); // return new fit as fitnum[0]=fitnum[0]+1
              if(finalden>minden){  // return the overall best rotation and translation matrixs
                   iass =2;  //only look for unassigned fits
                   findoneidentical(&iquit,del,iass,fitnum[0],ident);
                   for(ii=1;ii<ident[0][0][0];i++) // assign ident fit va1l = 10000 so not assigne to segment 1
                      fit[ident[ii][0][0] ].val1 = 100000;
              }
              fitnum[0]=origfitnum; // remove last addition
           }
           else{ // no more fits that can be assigned
             goto quitwhile2;
           }
        } // not assigned   
      }//i  
    }// while 
    quitwhile2:;
#endif        
   if(verbose>=1){
      printf(" ASSIGNED following fits to seg. # %d:",count);
      for(i=1;i<=nseg[0];i++)
         printf(" %d",nseg[i]);
      printf("\n");
   }
    iass=1; // 1jun01 - changed to this value test only fits witn sequence assignee in findoneidentical
//    iass=0; //  test all fits  in findoneidentical
int ifind,iplus,savsym[NSYM][2],ufit;
    ifind =1; // return best rot matrix after finding first sym mate
    iplus =0;
    for(ii=1;ii<=nseg[0];ii++){ 
       for(ufit=1;ufit<=origfitnum;ufit++){
          if( (fit[ufit].val1<=0) && (fit[ufit].val2> -900) ){
             symonefitass(&ifind,geo,den1,nseg[ii],savsym,&nsym,
                   brot,btran,del,ufit,&minden);  //ifind not used anymore
             if(minden>=minden0){ // if find sym operator - apply to all fits in seg1
               count = count+1;
               addone =0;
               for(jj=1;jj<=nseg[0];jj++){
                 refinerottran(geo,den1,brot,btran,maxside,rotm,tran,nseg[jj],
                      minden0,&finalden,&frout); // return new fit as fitnum[0]=fitnum[0]+1
//savefit();
                 if(finalden>minden0){  // return the overall best rotation and translation matrixs
                    findoneidentical(&iquit,del,iass,fitnum[0],ident);
                    if(verbose>=1){
                        printf("For afit=%d ufit=%d count=%d: found %d identical fits:",
                           nseg[jj],ufit,count,ident[0][0][0] );
                        for(j=1;j<=ident[0][0][0];j++)
                           printf(" %d",ident[j][0][0]);
                        printf("\n");
                   }
                   for(j=1;j<=ident[0][0][0];j++){
                       if(fit[ ident[j][0][0] ].val2 > -900){ // segment not assigned
                          addone = 1;
                          fit[ ident[j][0][0] ].val2 = -900 - count;
printf(" Assigning fit %d to segment %d\n",ident[j][0][0],count);
                       }
                    }
                  } // >minden
                  fitnum[0]=origfitnum; // remove last addition
                }//jj
                if(addone==0) // if did not another segment
                    count = count-1;
              } //minden
              else if( minden> -1000.0) { // Indicates fits  overlap
                 printf("ERROR: FIT %d  OVERLAPS %d - BUT POOR minden - FIT %d PROBABLY INCORRECTLY ASSIGNED\n",ufit,nseg[ii],ufit);
              }
          }// if seg not assigned
       }//ufit
    }//ii
#if 0
    for(isym = 1;isym <=nsym;isym++){  // Apply NCS fits in  seg = 1 = amax
       count = count+1;
       addone = 0;
       for(ii=1;ii<=nseg[0];ii++){ 
           amax = nseg[ii]; 
           refinerottran(geo,den1,brot[isym],btran[isym],maxside,rotm,tran,amax,
                minden,&finalden,&frout); // return new fit as fitnum[0]=fitnum[0]+1
           if(finalden>minden){  // return the overall best rotation and translation matrixs
               findoneidentical(&iquit,del,iass,fitnum[0],ident);
printf(" isym=%d amax=%d count=%d,finalden=%5.3f,frout=%5.3f iden[0]=%d,ident1=%d\n",
                    isym,amax,count,finalden,frout,ident[0][0][0],ident[1][0][0]);
savefit();
//exit(0);
               fitnum[0]=origfitnum; // remove last addition
               for(j=1;j<=ident[0][0][0];j++){
                 if(fit[ ident[j][0][0] ].val2 > -900){ // segment not assigned
                    addone = 1;
                    fit[ ident[j][0][0] ].val2 = -900 - count;
                 }
               }
           } // >minden
           fitnum[0]=origfitnum; // remove last addition
        } //ii
        if(addone==0) // if did not another segment
           count = count-1;
   }// isym
#endif
   for(i=1;i<=origfitnum;i++){ // By this point - all assigned fits should have segment assigned - if not 
      if( (fit[i].val1<0) && (fit[i].val2> -900) ){
         printf("WARNING:assigned fit %d  (val1=%d val2=%d) does not have seg. assigned\n",i,fit[i].val1,fit[i].val2);
         count=count+1;
         printf("   Assigning it segment # %d\n",count);  
         fit[i].val2 =  -900 - count;
      }
   }
   for(i=1;i<=origfitnum;i++){ // assign the remaining unassigned fits 
      if( (fit[i].val1>=0) && (fit[i].val2> -900) ){
         count=count+1;
         printf("   Assigning unassigned fit %d to segment # %d\n",i,count);  
         fit[i].val2 =  -900 - count;
      }
   }
   char segfile[300];
   strcpy(segfile,maidfile.fitfile);
   strcat(segfile,"_");
   strcat(segfile,callfile);
   strcat(segfile,"_seg");
   fp= fopen(segfile,"w");
   fprintf(fp,"%d\n",origfitnum); // total  # of fits
   for(i=1;i<=count;i++){
    maxres = -10;
    while(1){ // find lowest numbered res - use goto to exit
      minres = 100000;
      for(ii=1;ii<=origfitnum;ii++){  // save fits - grouped as segment, lowest residue first
         if( (fit[ii].val2 == (-900-i) )&&(fit[ii].val1<0) ){
            if( (-fit[ii].val1 <minres)&&( -fit[ii].val1>maxres) ){
               minres = -fit[ii].val1;
               imin=ii;
            }
         }
      }
      if(minres< 90000){//found a fit to save
         saveonefit(fp,fit,fitnum[imin],fitnum[imin]);
         maxres = minres;
      }
      else
         goto nexti2;   
    }// while
    nexti2:;
    for(ii=1;ii<=origfitnum;ii++)  // save the unasigned fits - grouped as segment
       if( (fit[ii].val2 == (-900-i) )&&(fit[ii].val1>=0) )
         saveonefit(fp,fit,fitnum[ii],fitnum[ii]);
  }//i
  fclose(fp);
  readfit(geo,segfile); // Need to read becuause above reorded fits
  strcpy(fitpdbfile,maidfile.fitfile); // Save the fit as fitfile.pdb
  strcat(fitpdbfile,"_");
  strcat(fitpdbfile,callfile);
  strcat(fitpdbfile,".pdb");
  writefittopdbgen(den1,geo,fit,fitnum,fitpdbfile);
  quitrout:;
  printf("FINISHED ROUTINE assignseqment()\n");
}
#endif

void rotatepdb(struct pdbfile& pdb,char pdbnew[])
/* Routine to rotate the screen pdb  to file pdbnew.
USES the first rottran operator in name.maidsymrottran file
*/
{
   int nsym;
   float rot[NSYM][3][3],tran[NSYM][3],xnew[3];
   FILE *fp,*fpout;
   char lresnumch[NCHAR]; 
   char dataline[200],atomch[30],atomnumch[NCHAR],atomname[NCHAR];
   char resname[NCHAR],seg[NCHAR];
   float sigma,bfact,x[3];

  printf("APPLYING ROT/TRAN OPERATOR #1 TO PDB - OUTPUT %s\n",pdbnew); 
  if( (fp=fopen(pdb.name,"r") ) ==NULL) {
      fprintf(stderr,"error opening file %s\n",pdb.name);
      exit(0);
  }
  fpout = fopen(pdbnew,"w");
  readnsym(maidfile.nsym,&nsym,rot,tran); // read rot/tran operators
  while(  fgets(dataline,199,fp) != NULL) {
    if(strncmp("ATOM",dataline,4)!=0){ // If does not begin with ATOM, just write again
       fprintf(fpout,"%s",dataline); 
    }
    else { // begins with ATOM
         if(dataline[21]!=' ') /* then standard pdb format */
            sscanf(dataline,"%s%s%s%s%s%s%8f%8f%8f",atomch,atomnumch,atomname,
               resname,seg,lresnumch, &x[0],&x[1],&x[2]);
        else{ /* xplor format */
           sscanf(dataline,"%s%s%s%s%s%8f%8f%8f%6f%6f%s",atomch,atomnumch,atomname,resname,lresnumch,
                                     &x[0],&x[1],&x[2],&sigma,&bfact,seg);
        }
        rottranvec(rot[1],tran[1],x,xnew);// Use first listed operator
        fprintf(fpout,"ATOM  %5s  %-3s %3s  %4s     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
                     atomnumch,atomname,resname,lresnumch,
                     xnew[0],xnew[1],xnew[2],sigma,bfact,seg);
        }
     }    
   fclose(fp);
   fclose(fpout);
}
