/* extendfit.c++ -15FEB2000  STARTED COMPLETELY NEW FILE -
    USING NEW ROUTINES THAT ASSUME A.A. ASSIGNMENT IS KNOWN.
  New Routines:

  I) For each fit, finds edited bones, and best assignment value,while
      varying gap by +/- 1.  If best value is much better than second 
      best, then assign sequence and start extending.
  II) After a.a assignment, assign residues and POWELL refine.
     Then for end 6 residues, run torsional dynamics to get best fit before
     extending.
  II) If next residue is not glycine, add one residue,trying all possible
         constraints - for each constraint - run torsional dynamics
      Find the best fit - and check distance from edited main.
      If it is glycine, then, try complete range of phi/psi and add another   
         residue, and continue until non-glycine.
           
Type 1 = helix, 2 = sheet, 3  = left helix 4= epsilon 5 = delta  6= gammar 7 = gammal 
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/

#include "drawwindow.h"
#include "pdbtolev.h"
#include "editbones.h"
#include "math.h"
#include "trace.h" // includes a number of "extern .... " lines
#include "peptide.h"
#include "tordyn.h"
#include "extendfit.h"
#include "extend2.h"
#include "initialize.h"
#include "bone.h"
#include "rankside.h"
#include "assignseq.h"
#include "symmate.h"
#include "fittopdb.h"
#include "combine.h"

#include "string.h"
#include "stdlib.h"
#define srand48 srand

#define MAXFUSE 300 // max number of fits can keep track of fusing data
extern struct extbones ext; // All the data for edited bones used to extend one end of a given fit
extern struct maiddatafile maidfile;
extern struct tordyn tordata;
extern struct tordyn tordata2;
extern struct rotomer rot[30];// rotomer data for each amino acid
extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern int fitnum[NFIT];// see trace.c++ for definition
extern int verbose; //Defined in trace.c++  if = 1, print a lot of stuff
extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version
extern float fminden[10]; // defined in bone.c++ and  mapmenu.c++
extern char*  nameres[30];
int origfitnum;  //keeps track of of numer of fit in fit file
int nadd; // global variable  = # of residues to to add at each step
int fitbeg, fitend; // = beg and end values in copyregfittogeo,e tc.
int totadd;  // total # of res. added after reading original fit 
int startres;// Global variable - fit res # where start extenting
char extenddata[300];
struct pepfit fitlim; //global - becuase use in maxfitdistgeo() called from tordyn
int iextend; // for degugging
int iseed;  
FILE *fpext;

int dodebug=0;
int countfuse= 0; // global variable - counts fusions - writes set at time of fusion to fuseset

struct besttype  // stores phi psi for each addition of best den of each type
{
    float phi[MAXADD];
    float psi[MAXADD];
    float qual;
    int main[3];
    int nbad;
    int type[MAXADD];
};


void writeset(char setfile[],struct pepset& fitset);
void sidecadist(struct geometry& geo,float capos[3],int imain,int iside,
    int pt0,int nmax,int mpt[],float* mdist,float pos[3]);


void copygeotosav(struct geometry& geo,struct savegeo sav)
/* copies geo.pos and side and type data
*/
{
   int i;

   copypos(geo.iside,geo.nend[1]-geo.nend[0],geo.nend[0],geo.nend[0],
                   geo.pos,sav.spos);// save initial  geo.pos
   for(i=geo.nend[0];i<=geo.nend[1];i++) // copy type  data 
       sav.stype[i] = geo.rtype[i];
   for(i=geo.nend[0]+1;i<geo.nend[1];i++) // copy  side data 
       sav.sside[i] = geo.iside[i];
}
        
void copysavetogeo(struct geometry& geo,struct savegeo sav)
/* copies geo.pos and side and type data
*/
{
   int i;

   copypos(geo.iside,geo.nend[1]-geo.nend[0],geo.nend[0],geo.nend[0],
                sav.spos,geo.pos);// save initial  geo.pos
   for(i=geo.nend[0];i<=geo.nend[1];i++) // copy type  data 
       geo.rtype[i]=sav.stype[i];
   for(i=geo.nend[0]+1;i<geo.nend[1];i++) // copy  side data 
       geo.iside[i]= sav.sside[i];
}

void aveoncaden(struct griddata& den1,struct geometry& geo,
       int n0,float fpos[NAFIT][NAA][3],float *aveden)
//CALLED ONLY IN trace.c++ store ave density of o,n,Ca  
{
    *aveden=0.0;
    *aveden = *aveden+ aveatom(den1,geo,3,n0,fpos); // den of O
    *aveden = *aveden+ aveatom(den1,geo,4,n0+1,fpos); // den of N
    *aveden = *aveden+ aveatom(den1,geo,1,n0+1,fpos); // den of Ca
    *aveden = *aveden/3.0;
}

void avecocaden(struct griddata& den1,struct geometry& geo,
       int n0,float fpos[NAFIT][NAA][3],float *aveden)
//CALLED ONLY IN trace.c++ store ave density of C,O,and Ca  
{
    *aveden=0.0;
    *aveden = *aveden+ aveatom(den1,geo,3,n0,fpos); // den of O
    *aveden = *aveden+ aveatom(den1,geo,2,n0,fpos); // den of C
    *aveden = *aveden+ aveatom(den1,geo,1,n0,fpos); // den of Ca
    *aveden = *aveden/3.0;
}



void avedenreggrad(struct griddata& den1,struct geometry& geo,
       float pos[NAFIT][NAA][3],float *ave,float *avecb2,
       int *nbad,int rbeg,int rend,float *grad)
/*
 In addition to finding aveden, finds a weighed average of 
  gradient of density between neighboring  main chain atomsatoms on
*/
{
   int j,count,countcb2,iadd;
   float aden;
   float ca[MAXADD+2],c[MAXADD+2],n[MAXADD+2];

   *ave = 0.0;
   *avecb2 = 0.0;
   *nbad = 0;
   *grad = 0;
   geo.npoor = 0;
   count = 0;
   countcb2 = 0;
   iadd = rend - rbeg; // should = nadd
   if(rbeg < geo.nend[0])
     rbeg = geo.nend[0];
   if(rend > geo.nend[1])
     rend = geo.nend[1];
   for(j=rbeg;j<=rend;j++){  // Ca
     count = count +1;
     aden= aveatom(den1,geo,1,j,pos);
     *ave = *ave + aden;
     ca[j-rbeg]= aden;
     if(aden< geo.poorrigid)
        geo.npoor = geo.npoor+1;
     if(aden< geo.minrigid){
            *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,1);
        }
   }
   for(j=rbeg+1;j<=rend;j++){ // N
     count = count +1;
     aden = aveatom(den1,geo,4,j,pos);
     *ave = *ave + aden;
     n[j-rbeg]= aden;
     if(aden< geo.poorrigid)
        geo.npoor = geo.npoor+1;
     if(aden< geo.minrigid){
        *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,4);
     }
   }
   for(j=rbeg;j<=rend-1;j++){ 
     count = count +1;
     aden = aveatom(den1,geo,2,j,pos);// C
     *ave = *ave + aden;
     c[j-rbeg] = aden;
     if(aden< geo.poorrigid)
        geo.npoor = geo.npoor+1;
     if(aden< geo.minrigid){
        *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,2);
     }
     count = count +1;
     aden = aveatom(den1,geo,3,j,pos);// O
     *ave = *ave + aden;
     *avecb2 = *avecb2 + aden;  // avecb2 = aveden for O
     countcb2 = countcb2+1;
     if(aden< geo.minrigid){
        *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,3);
     }
   }// j
   *ave = *ave/count;
   *avecb2 = *avecb2/countcb2;
   for(j=0;j<iadd;j++) //find grad of den
     *grad = *grad+fabs(ca[j]-c[j])+fabs(c[j]-n[j+1])+fabs(n[j+1]-ca[j+1]);
   *grad = *grad/iadd;
}


void avedenregnadd(struct griddata& den1,struct geometry& geo,
       float pos[NAFIT][NAA][3],float *ave,float *avecb2,
       int *nbad,int rbeg,int rend,float nave[],int ocount)
/*
New if ocount = 1, then include O in nbad and den
   else do not
New returns nave[0,1,iadd] for each of the residues
  where iadd is defineded by rbeg, rend;
  Still returns ave (all residues) and nbad (all residues)
Count both Ca ends as part of residue: Ca-C-0-N-Ca
*/
{
   int j,count,countcb2;
   int iadd; 
   float aden,caave[2*MAXADD+2];

   *ave = 0.0;
   *avecb2 = 0.0;
   *nbad = 0;
   geo.npoor = 0;
   count = 0;
   countcb2 = 0;
   if(rbeg < geo.nend[0])
     rbeg = geo.nend[0];
   if(rend > geo.nend[1])
     rend = geo.nend[1];
   iadd = rend - rbeg; // should = nadd
   for(j=0;j<iadd;j++)
     nave[j] = 0.0;
   for(j=rbeg;j<=rend;j++){  // Ca
     count = count +1;
     aden= aveatom(den1,geo,1,j,pos);
     *ave = *ave + aden;
     caave[j-rbeg]= aden;
     if(aden< geo.poorrigid)
        geo.npoor = geo.npoor+1;
     if(aden< geo.minrigid){
            *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,1);
        }
   }
   for(j=rbeg+1;j<=rend;j++){ // N
     count = count +1;
     aden = aveatom(den1,geo,4,j,pos);
     *ave = *ave + aden;
     nave[j-rbeg-1]= nave[j-rbeg-1]+aden;
     if(aden< geo.poorrigid)
        geo.npoor = geo.npoor+1;
     if(aden< geo.minrigid){
        *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,4);
     }
   }
   for(j=rbeg;j<=rend-1;j++){ 
     count = count +1;
     aden = aveatom(den1,geo,2,j,pos);// C
     *ave = *ave + aden;
     nave[j-rbeg]= nave[j-rbeg]+aden;
     if(aden< geo.poorrigid)
        geo.npoor = geo.npoor+1;
     if(aden< geo.minrigid){
        *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,2);
     }
if(ocount==1){ // Include )
     count = count +1;
     aden = aveatom(den1,geo,3,j,pos);// O
     *ave = *ave + aden;
     nave[j-rbeg]= nave[j-rbeg]+aden;
     *avecb2 = *avecb2 + aden;  // avecb2 = aveden for O
     countcb2 = countcb2+1;
     if(aden< geo.poorrigid)
        geo.npoor = geo.npoor+1;
     if(aden< geo.minrigid){
        *nbad = *nbad+1;
//fprintf(fpout," nbad for res = %d  atom = %d\n",j-geo.nend[0]+1,3);
     }
} // ocount
   }// j
   *ave = *ave/count;
   *avecb2 = *avecb2/countcb2;
   for(j=0;j<iadd;j++) // add Ca density at both ends to nave
     nave[j]=nave[j]+caave[j]+caave[j+1];
   for(j=0;j<iadd;j++) 
      nave[j] = nave[j]/5;
#if 0
   for(j=cbbeg;j<=cbend;j++)
      printf("cbden[%d]=%d  iside = %d",j,geo.cbden[j],geo.iside[j]);
printf("\nrbeg = %d rend = %d geo0,1 =%d %d usecb = %d CAend[x] =%5.4f  lastden = %5.4f\n",
     rbeg,rend,geo.nend[0],geo.nend[1],geo.usecb,pos[geo.nend[1]][1][0],*ave);
#endif
}
        
void rotonca(struct geometry& geo,float pos[NAFIT][NAA][3],float fpos[NAFIT][NAA][3],
      int n0,float rotv[3],float orig[3],float ang) //rotate o,n,and ca about rotv
// Rotate C, O, N and Ca by ang
{
   float rom[3][3];

   rotaxis(geo.de,geo.ep,ang,rotv,rom); /* rom is rotation matrix for rotation
         of chi about unit vector u2 which is assumed to be at origin*/  
   rotatesetnew(orig,pos[n0][3],fpos[n0][3],rom); // rotate O
   rotatesetnew(orig,pos[n0+1][4],fpos[n0+1][4],rom);// rotate N
   rotatesetnew(orig,pos[n0+1][1],fpos[n0+1][1],rom);// rotate Ca
}
        
void nrotcoca(struct geometry& geo,float pos[NAFIT][NAA][3],float fpos[NAFIT][NAA][3],
      int n0,float rotv[3],float orig[3],float ang) //rotate o,Cand ca about rotv
//Nterm -  Rotate C, O, and Ca by ang
{
   float rom[3][3];

   rotaxis(geo.de,geo.ep,ang,rotv,rom); /* rom is rotation matrix for rotation
         of chi about unit vector u2 which is assumed to be at origin*/  
   rotatesetnew(orig,pos[n0-1][3],fpos[n0-1][3],rom);// rotate O
   rotatesetnew(orig,pos[n0-1][1],fpos[n0-1][1],rom);// rotate Ca
   rotatesetnew(orig,pos[n0-1][2],fpos[n0-1][2],rom);// rotate C
}
 
void rotncb(struct geometry& geo,float pos[NAFIT][NAA][3],float fpos[NAFIT][NAA][3],
     int n0,float rotv[3],float orig[3],float ang)
 //rotate c and cb by ang about rotv, (origin) 
{
   float rom[3][3];

   rotaxis(geo.de,geo.ep,ang,rotv,rom); /* rom is rotation matrix for rotation
         of chi about unit vector u2 which is assumed to be at origin*/  
   rotatesetnew(orig,pos[n0][5],fpos[n0][5],rom); // rotate Cb
   rotatesetnew(orig,pos[n0][4],fpos[n0][4],rom);// rotate N
}

void rotccb(struct geometry& geo,float pos[NAFIT][NAA][3],float fpos[NAFIT][NAA][3],
     int n0,float rotv[3],float orig[3],float ang)
 //rotate c and cb by ang about rotv, (origin) 
{
   float rom[3][3];

   rotaxis(geo.de,geo.ep,ang,rotv,rom); /* rom is rotation matrix for rotation
         of chi about unit vector u2 which is assumed to be at origin*/  
   rotatesetnew(orig,pos[n0][5],fpos[n0][5],rom); // rotate Cb
   rotatesetnew(orig,pos[n0][2],fpos[n0][2],rom);// rotate C
}

void saverankfit(char rankedfile[],int nrank,struct pepfit fits[])
{
   int i;
   FILE *fp;

fprintf(fpout,"Saving %d ranked fit files in %s  length = %d\n",
       nrank,rankedfile,fits[1].nres);
   fp= fopen(rankedfile,"w");
   fprintf(fp,"%d\n",nrank); // total  # of ranked fits
   if(nrank>0) // print number of previous rankings on next line
       fprintf(fp,"%d\n",fits[1].irank[0]);
   for(i=1;i<=nrank;i++)
      saveonefit(fp,fits,i,i);
   for(i=1;i<=nrank;i++) // print value of main for each fit  
      fprintf(fp,"  %d %d %d",fits[i].main[0],fits[i].main[1],fits[i].main[2]);
   fprintf(fp,"\n");
   fclose(fp);  
}

void readfittemp(struct geometry& geo)
/* REads a specific file - not the fitfile - so does not interfere with 
   programs running trace in batch mode.
 read each fit into consecutively numbered bines
 If fitnum[i] is not consecutive - MAKE it consecutive
*/
{
   int i,irank0;
   FILE *fp;
   char fitfile[300];
//   Widget widget; 
//strcpy(fitfile,"fit_good.2mar99");
//strcpy(fitfile,"fit_good.29apr99"); //NOTE - DOES NOT EXIST  - WRITTEN OVER
strcpy(fitfile,"fit_good.21may99");
if(verbose>=1)
  fprintf(fpout,"Reading fitfile %s\n",fitfile);
 if( (fp=fopen(fitfile,"r")) == NULL){
     fitnum[0] = 0; // indicates no fit residues
#if 0  // do not use error message because it requires screen functions - not use in batch mode
         errormessage(toplevel,"ERROR: There are no fit files");
#endif
     fprintf(fpout,"ERROR: There are no fit files\n" );
 }
 else{ //fitfile opened
   fscanf(fp,"%d",&fitnum[0]); // total  # of fits;
   for(i=1;i<=fitnum[0];i++){ // each fit
     irank0=0;  // there are no rankings for files read in readfit
     readonefit(fp,fit,irank0,i);
   }// each fit 
   fclose(fp);  
 } // fitfile exists 
}

void readset(char setfile[],struct pepset& fitset)
/* First reads setfile which described set of rankedfits,
  Then reads each rankedfit consequtively into fit[]
*/
{
   int i,ii,j,nfits,nrank,count;
   FILE *fp;
   char rankfile[10],sfit[10],ranknum[10];

fprintf(fpout,"Reading set %s of rankedfits\n",setfile);
   fp = fopen(setfile,"r");
   strcpy(ranknum,setfile);
   strcat(ranknum,"fit");
   fscanf(fp,"%d",&fitset.fit[0]); // Number of initial fit files
   fscanf(fp,"%d",&fitset.rank[0][0]);//Number of ranked fit files (should = # fit files)
//   nfits = 0;
   nfits = fitnum[0];
   for(i=1;i<=fitset.rank[0][0];i++){ // File name of it ranked st
      fscanf(fp,"%d",&fitset.fit[i]);
      strcpy(rankfile,ranknum);
      itoa(fitset.fit[i],sfit);
      strcat(rankfile,sfit); 
      readrankfit(rankfile,&nrank,fit,nfits+1); // means start with fits[nfit]            
      fitset.rank[i][0]=nrank; // # of fits in ith rankedfitg
      for(j=1;j<=nrank;j++) // # of fit file corresponding to ranked fits in ith set
         fitset.rank[i][j]=nfits+j;
      nfits = nfits+nrank;
      count=0;
      for(ii=nfits-nrank+1;ii<=nfits;ii++){// Assign rankfit data
         count = count+1;
         fit[ii].rankfit[0]=fitset.fit[i];
         fit[ii].rankfit[1]=count;
      }
   }
   fitnum[0]=nfits;
}

void writeset(char setfile[],struct pepset& fitset)
{
  int i;   
  FILE *fp;
  char ranknum[10];
   
fprintf(fpout,"Saving set file %s \n",setfile);
   fp = fopen(setfile,"w");
   strcpy(ranknum,setfile);
   strcat(ranknum,"fit");
   fprintf(fp,"%d\n",fitset.fit[0]); // Number of initial fit files
   fprintf(fp,"%d\n",fitset.rank[0][0]);//Number of ranked fit files (should = # fit files)
   for(i=1;i<=fitset.rank[0][0];i++)
     fprintf(fp," %d",fitset.fit[i]);
   fprintf(fp,"\n");
   fclose(fp);
}



#if 1 // new version  - calls linedist
void maincadist(struct geometry& geo,float capos[3],int imain,int pt0,int nmax,
    int mpt[],float* mdist,float pos[3],float *bdist) 
/*
 NEW VERSION IS NOT Recursive function -do not need to call again if all points on one main
 Starting with pt0, move out along mains for a total of nmax steps
   and find closest point
 NEW:  closepoint to be with the p1 to p2 interval, 
    NOTE: since try both intervals on each side of closepoint are
      tried, this allows either pt either branch to be used
   (i.e.  0<a< |v}) (see routine linedist()
  Allow some leaway;  -0.5 < a < |v|+0.5
linedist returns idist =0 if not within this interval
     = 1 if it is   
Try both lower and upper intervals
If a not in this interval for for either interval
   then return simple dist from pt to closest bonept.
 Function returns npt[2] and mdist
   mpt[0] = main num of close pt
   mpt[1] = close pt #
   mpt[2] = 0 if close pt on mainbone, = side # if close pt on side that links to this main
   mdist  = close dist
   pos[3] = postion of interpolated closest distance
NOTE:   count runs from 1,2,... nmax corresponding to pt0,pt0+1,pt0+nmax-1
NO  - NOT IN NEW VERSION ALSO - MUST ALLOW TO MOVE OUT ON SIDE CHAINS THAT CONNECT TO OTHER MAINS 
  if ctype[0]=1 for a point, then call sidecadist, very similar to main ca dist 
*/
{
   int i,k,count,minpt;
   float cdist[40]; //10 must be greater than nmax*nadd+1
   float mindist,tpos[3];
   int idist1,idist2;
   float mdist1,mdist2,mindist2;

   if(pt0<0)
     pt0=0;
   i= pt0;
   count = 0;
   mindist = 2000.0;
   cdist[0] = 1000.0; // set to large value, - 0 pt cannot be near best
   while( (i<ext.mainone.num)&&(count<nmax) ){
     count = count+1; // count runs from 1,2,... nmax
     cdist[count] = distv(capos,ext.mainone.mainpt[i].pos);
     if( cdist[count]<mindist){
        mindist = cdist[count];
        minpt = i;
//        cmin = count;
     }
     i = i+1;
   }  // while
   *bdist = 0.0;
   for(i=pt0+1;i<=minpt;i++){ // find distance along bones to minpt
      *bdist = *bdist +distv(ext.mainone.mainpt[i].pos,ext.mainone.mainpt[i-1].pos);
   }
#if 0 // New  4oct00 - now applied in calling routine
   if(*bdist > 7.5){// NEW 20sep00- prevents from jumping around loop
      *mdist = 1000.0; // give it a bad value
      goto quitrout;
   }
#endif
   mindist2 = 2000.0;
   if(minpt-1<0){
      mdist1 = 1000.0;
      idist1 = 0; // indicates not in interval
   }
   else
      linedist(ext.mainone.mainpt[minpt-1].pos,ext.mainone.mainpt[minpt].pos,
               capos,&mdist1,pos,&idist1);  // try lower interval 
   if(idist1 ==1) { // in interval
     mindist2 = mdist1; 
     for(k=0;k<=2;k++)
        tpos[k] = pos[k];
   }
   linedist(ext.mainone.mainpt[minpt].pos,ext.mainone.mainpt[minpt+1].pos,
           capos,&mdist2,pos,&idist2); // upper interval
   if( (idist2 ==1)&&(mdist2 < mindist2) ) { // then use lower interval
      *mdist = mdist2;
   }
   else if(idist1==1){ // use first interval
      *mdist = mdist1;
      for(k=0;k<=2;k++)
        pos[k] = tpos[k]; 
   }
   else{ // not in either interval
     *mdist = mindist; 
      for(k=0;k<=2;k++)
        pos[k] =ext.mainone.mainpt[minpt].pos[k]; 
   }
//printf("bonestep = %d mdist = %5.3f  bonedist = %5.3f\n",minpt-pt0,*mdist,bdist); 
   mpt[0]= imain;
   mpt[1] = minpt;
   mpt[2]= 0; // indicates not on connecting side
#if 0 // debug 
   float tdist = distv(capos,pos);
   if(fabs(tdist- *mdist)> 0.001)
     printf("ERROR:  tdist = %5.3f  *mdist = %5.3f\n",tdist,*mdist);
#endif
   quitrout:;
}
#endif


void copyfpostogeo(int nend,float fpos[NAFIT][NAA][3],struct geometry& geo)
{
  int k;           
           for(k=0;k<=2;k++){
              geo.pos[nend][3][k] = fpos[nend][3][k];
              geo.pos[nend+1][4][k] = fpos[nend+1][4][k];
              geo.pos[nend+1][1][k] = fpos[nend+1][1][k];
           }
}

#if 1
void checkoverlap(struct geometry& geo,struct pepfit& fit1,
               struct pepfit& fit2,int *nlap,float *maxdist,float *mindist)
/*
  NEW 8dec99:    Have completely rewritten using same procedure as in linedist().
      Now: finds value of a (dist along Ca-Ca) and perpdist for each
        Ca.  Find best perpdist and calculates overlap.
  Checks if fit1 and fit2 overlap enough to try and fuse
    and then return the  maxdist apart for all the  atoms that art to be fused
  Define fit 2 as original fit that is met
    and fit1 as the extendfit that want to fuse with fit2
  Returns nlap:  = -1 if no overlap,
               and = j  if j residues overlap, that is, want to remove
               j residues from fit1 before trying to fuse.
  Define: cend[] = pos of nterm Ca in fit2.
          c2[1 .. 5][] = pos of cterm Ca of fit1(=1) and next for Ca' of fit2
  Then: For each unit vector u = c2[i+1]-c2[i], find the mag a defined
    such that the perpendicular projection (=v)  to u meets u at postion 
        c2[i]+ au.
Find as as follows:  
      v = cend - (c2[i]+ a u)  and since u dot v = 0
      a = (cend dot u) - (c2[i] dot u)
The if   cadist/2 < a < 1.5 cadist, then 
    fit overlaps the c2[i] to c2[i+1] residue
   and should return nlap = i-1;
  For example, if overlaps c2[1] to c2[2], then nlap = 0, and do not want remove any
     residues for fusing.
NOTE:  After extension, new fits are stored in fit[origfitnum+i]
   and fit it meets = ext.metend
*/
{
  int i,j,k,iend;
  float a,mag,u[3],cend[3],c2[8][3],dist;
  float p[3],magp,perp,minperp,mina;
  int mini;

  *nlap = -1;  // idicates no overlap
  *maxdist =  -1.0;
  *mindist  = 1000.0;
   mina = -1000.0;
//New routineusing simpler routine -see void linedist() for details
  for(k=0;k<=2;k++){// x,y,z
    cend[k] = fit2.pos[0][1][k]; // nterm end Ca of fit2
    for(j=1;j<=6;j++) // C term Ca atoms in fit1
      c2[j][k] = fit1.pos[fit1.nres-j][1][k];
  } // k
  minperp = 1000.0;
  for(i=1;i<=5;i++){ // 4 intervals
    for(k=0;k<=2;k++)
      u[k] = c2[i][k]-c2[i+1][k];
    mag = sqrt( dot(u,u) );
    for(k=0;k<=2;k++)
      u[k] = u[k]/mag;
    for(k=0;k<=2;k++)
      p[k]= cend[k] - c2[i+1][k];
    magp = sqrt(dot(p,p));
    a= dot(u,p);
    perp = sqrt((magp*magp - a*a) );
    if( (a> -0.5*geo.cadist)&&(a<= geo.cadist+0.5*geo.cadist) ){
       if(perp<minperp){
         minperp = perp;
         mini = i;
         mina = a;
       }
     }
//printf("i = %d, a = %5.3f  perpdist = %5.3f\n",i,a,perp);
  }
  if( (mina> -0.5*geo.cadist)&&(mina<= 0.5*geo.cadist) )
       *nlap = mini-1;
  else if( (mina> 0.5*geo.cadist)&&(mina<= geo.cadist+0.5*geo.cadist) )
          *nlap = mini-2;
//  printf("nlap = %d\n",*nlap);

  quitrout:;
  if( *nlap>=0){ // then find maxdist after remove nlap res. from fit 
    iend = fit1.nres -1- *nlap;  // end res number of fit1
    dist = distv(fit1.pos[iend][1],fit2.pos[1][1]); // next to last Ca
    if(dist> *maxdist)
       *maxdist = dist;
    if(dist< *mindist)
       *mindist = dist;
    dist = distv(fit1.pos[iend-1][1],fit2.pos[0][1]); //end Ca
    if(dist> *maxdist)
       *maxdist = dist;
    if(dist< *mindist)
       *mindist = dist;
    dist = distv(fit1.pos[iend][4],fit2.pos[1][4]); // N 
    if(dist> *maxdist)
       *maxdist = dist;
    if(dist< *mindist)
       *mindist = dist;
    for(i=2;i<=3;i++){ // C, O 
        dist = distv(fit1.pos[iend-1][i],fit2.pos[0][i]); 
        if(dist> *maxdist)
          *maxdist = dist;
        if(dist< *mindist)
          *mindist = dist;
    }
  } // if nlap
if(verbose>=1)
    printf("In Checkoverlap nlap = %d,maxdist = %5.3f mindist = %5.3f\n",*nlap,*maxdist,*mindist);
}   
#endif

void fuseextfit(struct griddata& den1,struct geometry& geo,int nlap,
    struct pepfit fitext,struct pepfit fitmet,int metend)
/*
New version - copies only regions of fit to geo
 1)Copies from istart to end of fitext to geo
 and from 0 to some number of fitmet
 2) After dynamics
     Just copy fitmet to the end of fitext
 3) then use copyreggeo to copy geo to the middle of this extend fitmet

 Fuse fit-ext (the new extended fit - cterm end) to fit-met 
    after removing nlap residues from fitext
  When saving, delete the lowered numbered fit
*/ 
{
      FILE *fp;
      int i,icount,inum,i1,i2;
      int istart,ntotal;
      struct pepfit fittemp;
      int origendmove;
      float origrfor,origfuseconst;

#if 0 // for debugging - just save extendfit and quit
   copyfittofit(fitext,fit[geo.fitn],0,fitext.nres);
   savefit();
   exit(0);
#endif
      origrfor = tordata.rfor;
      origfuseconst = tordata.fuseconst;
      tordata.rfor=0.0; // 16mar01 - Do not use random force during fusion
      origendmove = geo.endmove;
      geo.fuse[0]=2;
      geo.fixend = 3;  // Indicates dynamcis about center
      geo.quit = 0;
      fitext.nres = fitext.nres-nlap;  // Remove nlap res for Cterm end
      copyreg2fitstogeo(geo,fitmet,fitext,&istart); // this also adds fitmet to end of fitext
//      copy2fitstogeo(geo,fitmet,fitext);
      if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
         tordata.adata[1] = 2; 
      if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
         tordata2.adata[1] = 2;

//geo.endmove = 4;

#if 1 //New version - slowly increase endmove and increase fuseconst
      tordata.fuseconst = origfuseconst/3.0; //16mar01 - run twice - first with low fuseconst
      geo.endmove = 1;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("First fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = 3;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Second fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = origendmove;
      for(i=1;i<=30;i++){ // run tordyn 20 times
         if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Third fuseconst = %5.3f: Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      tordata.fuseconst = origfuseconst; //16mar01 - run twice - first with low fuseconst
      geo.endmove = 1;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("First fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = 3;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Second fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = origendmove;
      for(i=1;i<=30;i++){ // run tordyn 20 times
         if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Third fuseconst = %5.3f: Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
#endif

      endtordyn:;
      copyavefpostopos(geo,geo.fnum);

      ntotal = fitext.nres + fitmet.nres;
      fittemp.nalloc = -10;
      allocatefit(fittemp,ntotal+2); // make sure big enough allocation
      addfittofit(fitext,fitmet,fittemp); // this assigns fittemp.val = fitext.val
      if(fitmet.val2== -3){//SPECIAL CASE: extended fit met a fit that could not be extended in nterm direc.
            if(verbose>=1)
               printf("SPECIAL CASE:extended fit met a fit that could not be extended in nterm dir set val2 = -1\n");
            fittemp.val2= -1; // Means cannot extend in cterm dir. 
      }
      else
           fittemp.val2 = fitmet.val2; //if val2 = -1, cannot extend.
      fittemp.nres  = ntotal -2; // new value after addition of fit2
      copyreggeotofit(geo,fittemp,istart);

      fitnum[0] = origfitnum;
      if( geo.fitn < metend){
         i1 = geo.fitn;
         i2 = metend;
      }
      else{
         i1 = metend;
         i2 = geo.fitn;
      }
      fp= fopen(maidfile.fitfile,"w");
      icount = 0;
      fprintf(fp,"%d\n",fitnum[0]-1); // total  # of fits
      for(i=1;i<=fitnum[0];i++){ //save fit, deleting lower number of the 2 fused fits and relacing higher number wit the fused fit
         if(i==i2){ // higher number
           icount = icount+1; 
           inum = fitnum[0]+1;
//           copygeotofit(geo,fit[inum],val); // Copy 
           copyfittofit(fittemp,fit[inum],0,fittemp.nres);
           fit[inum].irank[0]= 0; // This makes sure saveonefit will not write ranking
           saveonefit(fp,fit,inum,icount);
         } 
         else if (i==i1){ // skip lower number, do nothing
         }
         else {
           icount = icount+1;
           saveonefit(fp,fit,i,icount);
         }
      }
      fclose(fp);
      geo.fixend = 0;  // if fixed = 3, draw connect geo
      readfit(geo,maidfile.fitfile);
      freepepfit(fittemp);
  quitrout:;
  tordata.rfor= origrfor;
  tordata.fuseconst= origfuseconst;
  geo.endmove=origendmove;
}


void fuseextfit2(struct griddata& den1,struct geometry& geo,int nlap,
    struct pepfit fitext,struct pepfit fitmet,int metend)
/*
Called from combine.c++
Only difference is that fused fit is stored as higher numbered fit
  and lower numbered fit is NOT deleted
ALSO - IF geo.quit = 1, then no not fuse
   Happens if extend makes a new fit that is out of density.
*/ 
{
      int i,i2;
      int istart,ntotal;
      struct pepfit fittemp;
      int origendmove;
      float origrfor,origfuseconst;

      origrfor = tordata.rfor;
      origfuseconst = tordata.fuseconst;
      origendmove = geo.endmove;
      tordata.rfor=0.0; // 16mar01 - Do not use random force during fusion
      if(verbose>=1)
         printf("Fusing fits %d and %d\n",geo.fitn,metend);
      geo.fuse[0]=2;
      geo.fixend = 3;  // Indicates dynamcis about center
      geo.quit = 0;
      fitext.nres = fitext.nres-nlap;  // Remove nlap res for Cterm end
      copyreg2fitstogeo(geo,fitmet,fitext,&istart); // this also adds fitmet to end of fitext
//      copy2fitstogeo(geo,fitmet,fitext);
      if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
         tordata.adata[1] = 2; 
      if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
         tordata2.adata[1] = 2; 
#if 0
      for(i=1;i<=50;i++){ // run tordyn 50 times
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Quitting tordyn after %d runs because geo.quit =1\n",i); 
            printf("Quitting fuseextfit2 without fusing\n");
            goto quitrout;
         }
      }
#endif

#if 1 //New version - slowly increase endmove and increase fuseconst
      tordata.fuseconst = origfuseconst/3.0; //16mar01 - run twice - first with low fuseconst
      geo.endmove = 1;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("First fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            printf("Quitting fuseextfit2 without fusing\n");
            goto quitrout;
         }
       }//i
      geo.endmove = 3;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Second fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            printf("Quitting fuseextfit2 without fusing\n");
            goto quitrout;
         }
       }//i
      geo.endmove = origendmove;
      for(i=1;i<=30;i++){ // run tordyn 20 times
         if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Third fuseconst = %5.3f: Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
             printf("Quitting fuseextfit2 without fusing\n");
            goto quitrout;
         }
       }//i
      tordata.fuseconst = origfuseconst; //16mar01 - run twice - first with low fuseconst
      geo.endmove = 1;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("First fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            printf("Quitting fuseextfit2 without fusing\n");
            goto quitrout;
         }
       }//i
      geo.endmove = 3;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Second fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            printf("Quitting fuseextfit2 without fusing\n");
            goto quitrout;
         }
       }//i
      geo.endmove = origendmove;
      for(i=1;i<=30;i++){ // run tordyn 20 times
         if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Third fuseconst = %5.3f: Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            printf("Quitting fuseextfit2 without fusing\n");
            goto quitrout;
         }
       }//i
#endif
      copyavefpostopos(geo,geo.fnum);

      ntotal = fitext.nres + fitmet.nres;
      fittemp.nalloc = -10;
      allocatefit(fittemp,ntotal+2); // make sure big enough allocation
      addfittofit(fitext,fitmet,fittemp); // this assigns fittemp.val = fitext.val
      fittemp.val2 = fitmet.val2; //if val2 = -1, cannot extend.
      fittemp.nres  = ntotal -2; // new value after addition of fit2
      copyreggeotofit(geo,fittemp,istart);
      if( geo.fitn < metend){
         i2 = metend;
      }
      else{
         i2 = geo.fitn;
      }
      if(i2>=NFIT){
        if(verbose>=1)
          printf("ERROR:  i2 (=%d) >=NFIT in fuseextfit2- quitting and not using fused\n",i2);
        geo.quit=2; // Indication quit without fusing nad reached limited of NFIT
        goto quitrout2; // Need to freepepfit
      }
      if(verbose>=1)
         printf("Copying fused fit to fit %d\n",i2);
      copyfittofit(fittemp,fit[i2],0,fittemp.nres); // store as higher numbered fit
      quitrout2:;
      freepepfit(fittemp);
  quitrout:;
  tordata.rfor= origrfor;
  tordata.fuseconst= origfuseconst;
   geo.endmove=origendmove;
}


void fuseextfit3(struct griddata& den1,struct geometry& geo,
    struct pepfit fitext,struct pepfit fitmet,struct pepfit& newfit)
/*
Called from trace - after geo fit that is being extended met another fit
  fuse two fits - and returns the fused fit newfit
*/ 
{
      int i;
      int istart,ntotal;
      struct pepfit fittemp;
      int origendmove;
      float origrfor,origfuseconst;

      origrfor = tordata.rfor;
      origfuseconst = tordata.fuseconst;
      tordata.rfor=0.0; // 16mar01 - Do not use random force during fusion
      origendmove = geo.endmove;
      geo.fuse[0]=2;
      geo.fixend = 3;  // Indicates dynamcis about center
      geo.quit = 0;
      copyreg2fitstogeo(geo,fitmet,fitext,&istart); // this also adds fitmet to end of fitext
//      copy2fitstogeo(geo,fitmet,fitext);
      if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
         tordata.adata[1] = 2; 
      if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
         tordata2.adata[1] = 2; 
#if 0
      for(i=1;i<=50;i++){ // run tordyn 30 times
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Quitting tordyn after %d runs because geo.quit =1\n",i); 
            goto endtordyn;
         }
      }
#endif
#if 1 //New version - slowly increase endmove and increase fuseconst
      tordata.fuseconst = origfuseconst/3.0; //16mar01 - run twice - first with low fuseconst
      geo.endmove = 1;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("First fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = 3;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Second fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = origendmove;
      for(i=1;i<=30;i++){ // run tordyn 20 times
         if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Third fuseconst = %5.3f: Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      tordata.fuseconst = origfuseconst; //16mar01 - run twice - first with low fuseconst
      geo.endmove = 1;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("First fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = 3;
      for(i=1;i<=15;i++){ // run tordyn 15 times
         if(tordata.adata[1]>2) //NOTE: ESSENTIAL TO RESTART IF INCREASE  geo.endmove This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Second fuseconst = %5.3f:Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
      geo.endmove = origendmove;
      for(i=1;i<=30;i++){ // run tordyn 20 times
         if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
         if(tordata2.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata2.adata[1] = 2; 
         tordyn(den1,geo,tordata);
         if(geo.quit==1){
            printf("Third fuseconst = %5.3f: Quitting tordyn after %d runs because geo.quit =1\n",tordata.fuseconst,i); 
            goto endtordyn;
         }
       }//i
#endif
      endtordyn:;
      copyavefpostopos(geo,geo.fnum);
      ntotal = fitext.nres + fitmet.nres;
      fittemp.nalloc = -10;
      allocatefit(fittemp,ntotal+2); // make sure big enough allocation
      addfittofit(fitext,fitmet,fittemp); // this assigns fittemp.val = fitext.val
      if(fitmet.val2== -3){//SPECIAL CASE: extended fit met a fit that could not be extended in nterm direc.
            if(verbose>=1)
               printf("SPECIAL CASE:extended fit met a fit that could not be extended in nterm dir set val2 = -1\n");
            fittemp.val2= -1; // Means cannot extend in cterm dir. 
      }
      else
           fittemp.val2 = fitmet.val2; //if val2 = -1, cannot extend.
      fittemp.nres  = ntotal -2; // new value after addition of fit2
      copyreggeotofit(geo,fittemp,istart);
      copyfittofit(fittemp,newfit,0,fittemp.nres);
      freepepfit(fittemp);
  quitrout:;
  tordata.rfor= origrfor;
  tordata.fuseconst= origfuseconst;
   geo.endmove=origendmove;
}

void maxfitdist(struct geometry& geo,struct pepfit& fit1,int nlap,
        int *nlim)
/* Returns the number (nmax) of N and Ca that are more then limitdist from 
     edited main
 Only 1 edited main so use "1" in maincadist
*/
{
   int pt0,pt[3],ires;
   float ndist,cadist,pos[3];
   float limitdist = 1.0;
   float bdist;

   pt0 = 1; // main pt number where start checking for dist
   *nlim = 0;
   readeditbones("origedit"); // read the original set of bones
   for(ires = startres;ires < fit1.nres-nlap-1;ires++){ // do not count last residue since i is counted in maxdist
        maincadist(geo,fit1.pos[ires][4],1,pt0,geo.gnmax,pt,
                                 &ndist,pos,&bdist); // N dist
//printf(" For ires = %d ndist = %5.3f\n",ires,ndist);
        if(ndist > limitdist)
           *nlim = *nlim+1;
        maincadist(geo,fit1.pos[ires][1],1,pt0,geo.gnmax,pt,
                                 &cadist,pos,&bdist); // Ca dist
//printf(" For ires = %d   cadist = %5.3f\n",ires,cadist);
        if(cadist > limitdist)
           *nlim = *nlim+1;
        pt0 = pt[1];
   }
}


void maxfitdistgeo(struct geometry& geo,int *nlim)
/*
NOTE - NOT USED ANYMORE -
Called from tordyn:  First copies geo to fit, 
  then calls maxfitdist
*/
{
   copyreggeotofit(geo,fitlim,fitbeg);
   fitlim.nres = fitend+nadd+1;
//printf("geo.nend[1]=%d  geo.nend[0]=%d fitlim.pos[fitlim.nres-1][1][0] = %5.3f\n",geo.nend[1],geo.nend[0],fitlim.pos[fitlim.nres-1][1][0]);
   maxfitdist(geo,fitlim,0,nlim);// "0" = nlap so find dist to end of fit
//stcount = stcount+1;
//printf("stcount = %d nlim = %d fitbeg = %d  fitend = %d\n",stcount,*nlim,fitbeg,fitlim.nres);
//   freepepfit(fitlim);
}


void makeeditbones(struct griddata& den1,
       struct geometry& geo,struct bonedata& bone1,int ifit)
/*
  Makes edited bones at Cterm end (geo.fixend==2) or 
     nterm end (geo.fixend==1) of fit[ifit]
  This is the editbones part of old extendfit() routine.
       HOWERVER, BEWARE:
           The values of geo.main[0,1,2] are not NOW saved in fitfile. Thus
               if have previous extended fit and saved geo, and then reread
               fit and try to extend, will have WRONG geo.main[] values.
           Also - if save geo - same problem - cannot continue extending
*/ 
{
  int n0;

  ext.numlim = 40; // NOT USED NOW Maximum distance to trace main chain from origin
  ext.distlim = 23.0; // Max distance (Angstroms) trace chain from start
//  ext.distlim = 15.0; // Max distance (Angstroms) trace chain from start
  ext.sidelim=7; // If number of side pts > sidelim, then assume it is a main branch
  if(verbose>=2)
       fprintf(fpout,"Starting with previous fit # %d:  fit.nres = %d\n",ifit,fit[ifit].nres);
   ext.ifit = ifit;
   n0=fit[ifit].nres-1; // C term residue  N- Ca at end
   fit[ifit].irank[0]=0; // No previous rankings since not a rankedfit
   geo.main[0]=0; // Reset so that will start new extension.
   if(geo.fixend==2){ // extending at cterm end
//     editbones(den1,fit[ifit].pos[n0-1][2],fit[ifit].pos[n0][1],geo,bone1); // C . . Ca postion
       editbones(den1,fit[ifit].pos[n0-1][1],fit[ifit].pos[n0][1],geo,bone1); // Ca . . Ca postion
#if 0 // Now done in editbones for extendmain on fit
       if(ext.metend>0)// New  11jul00 - set first pt = Cterm Ca
          for(k=0;k<=2;k++) 
              ext.mainone.mainpt[0].pos[k] = fit[ifit].pos[n0][1][k];
#endif
   }
   else{ // extending at nterm
       editbones(den1,fit[ifit].pos[1][4],fit[ifit].pos[0][1],geo,bone1); // N . . Ca postion
//exit(0);
#if 0
       if(ext.metend>0)
         for(k=0;k<=2;k++) // New - set first pt = Nterm Ca
            ext.mainone.mainpt[0].pos[k] = fit[ifit].pos[0][1][k];
#endif
   }
   quitrout:;
}   


void batchextendcomb(struct griddata& den1,
          struct geometry& geo,int icon,int *foundcon)
/*
Modified version of batchextend, called from conbine.c++ -
  Looks for fits that are connected to other fit = icon,
 Returns foundcon = fit that it met,
   If no connection found - returns foudncon = -1
geo.fixend = 2  - cterm extension   =1  nterm ext.
*/
{
  int i,j,k,ires,setside;
  struct bonedata bone;  // this is the low density cutoff (0.8) bones
  int iden;
  int ibestden,connum,cnum,fitext,fitmet,maxngap;
  char ch[10];
  FILE *fp;

  *foundcon = -1;  // default, did not find connection
  geo.fitn = icon; 
  if(geo.nowset == 0){ // set not already allocated
       imat3(geo.set,den1.griddim[0],den1.griddim[1],den1.griddim[2]);// allocates array 3 dim array (mat)
       geo.nowset = 1;
  }
  for(iden = 6;iden>=0;iden--){// start with den[2.0] and reduce
    for(i=0;i<den1.griddim[0];i++){ // set all grid points =  0
     for(j=0;j<den1.griddim[1];j++)
       for(k=0;k<den1.griddim[2];k++)
          geo.set[i][j][k] = 0;  // set = 0
    }
    for(i=1;i<=fitnum[0];i++){ // set points for each fit - set side and end separately
        if(verbose>=2)
            fprintf(fpout,"Setting set[][][] for fit %d\n",i);
         geo.setval = geo.setside; // setside = setmain = -1

             setside = 2; // set 3 side pts
         for(ires=1;ires<fit[i].nres-1;ires++){ // New 29jun00
           setonesidefit(ires,den1,geo,fit[i],geo.set,setside);
         }// ires
         geo.setval = geo.setmain;
         setptsfit(den1,geo,fit[i],geo.set,0); // set points for main only - not side
         geo.setval = i;
         setendptsfit(den1,geo,fit[i],geo.set); // set points for end of main only - not side
    }//i
    bone.minden = fminden[iden]; 
    strcpy(bone.bonename,den1.mapname);
    strcat(bone.bonename,".bone"); // bone  name is mapname with .bone_minden appended
    sprintf(ch,"%1.1f",bone.minden);
    strcat(bone.bonename,ch);
    fp=fopen(bone.bonename,"r");
    if(fp == NULL){ // file does not exist - NOTE-closing NULL file cause crash in Linux
//         fclose(fp);
         findbone(den1,bone);
    }
    else{
       fclose(fp);
       readbonefile(bone); // This reopens bone file
    }
    makeeditbones(den1,geo,bone,geo.fitn);
    if(ext.metend>0){
       if( (fit[geo.fitn].val1<0)&&(fit[ext.metend].val1<0) ){ //both fitext and fitmet has sequence assignment - allow larger gap
         ibestden = 2; // this is passed to lowerconnectumca but not used in this routine
         lowerconnectnumca(geo,ibestden,&cnum,&maxngap); // should work for both cterm and nterm extension         
         if(geo.fixend==2){  // extend at cterm end
            fitext = geo.fitn;
            fitmet = ext.metend;
         }
         else{// extend at nterm end
            fitext = ext.metend;
            fitmet = geo.fitn;
         }
         printf("    fitext (%d)= %d and fitmet (%d)=%d\n",
                    fitext,-fit[fitext].val1,fitmet,-fit[fitmet].val1);
         connum = -fit[fitmet].val1 +fit[fitext].val1 - fit[fitext].nres +1;
         if( (connum>=cnum-4) && (connum<=cnum+4) ){ //Allow error of +/- 4 connum consistent with cnum1
            printf(" Gap number  = %d for assigment is consistent with bones\n",
                              connum);
            *foundcon = ext.metend;
         }
         else {
            printf(" Gap number  = %d for assigment is NOT consistent with bones\n",
                       connum);
         }
         *foundcon = ext.metend;
         goto quitrout;
      }// if both assigned
    } // ext.metend>0
  }// iden
  quitrout:;
}

void batchextend(struct griddata& den1,
          struct geometry& geo,int *foundedit,int *ibestden)
/*General Routine to extend fit
 NEW if cannot find edited bones, returns foundedit = 0
 NEW: only check one density value lower that the first that metend
 I) Finds the highest density that give extension that meets another fit
    1) Starts with  bone at den 1.6, make edit bones -
      If metend = 0, repeat at den 1.4
      Continue to minden = (say 1.0) - fit still metend =0 then quit
      NEW - use all denlevels, use the shortest main
    2) Using this best den that meets end, extend as far as possible
fminden[0,1,....] ={0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6}
*/
{
  int i,j,k,ires,setside;
  struct bonedata bone;  // this is the low density cutoff (0.8) bones
  struct mainonebone  mainonesav;
  int iden,minlength,testmetend,firstiden,savlength,savmetend;
  float bestden;
  char ch[10];
  FILE *fp;

  minlength = 10000;
  firstiden = -10;
  geo.userank = 0;
  *ibestden = -1; // default - no success
  *foundedit = 1; // default - found edited bones
  testmetend = 0; //This means when saveeditbones, to not filter
  if(verbose>=1)
     fprintf(fpout,"Starting with previous fit # %d:  fit.nres = %d\n",geo.fitn,fit[geo.fitn].nres);
//verbose=2;
//geo.miniden = 0;
  for(iden = 6;iden>=geo.miniden;iden--){//13feb01 - start with den[2.0] and reduce to geo.miniden (first call miniden = 3 - bone1.4)
//  for(iden = 6;iden>=6;iden--){//  just use 1 den
   for(i=0;i<den1.griddim[0];i++){ // set all grid points =  0
     for(j=0;j<den1.griddim[1];j++)
       for(k=0;k<den1.griddim[2];k++)
//          geo.set[i][j][k] = -10;
          geo.set[i][j][k] = 0;  // set = 0
   }
   for(i=1;i<=fitnum[0];i++){ // set points for each fit - set side and end separately
//        copyfittogeo(geo,fit[i]);       
        if(verbose>=2)
            fprintf(fpout,"Setting set[][][] for fit %d\n",i);
         geo.setval = geo.setside; // setside = setmain = -1

             setside = 2; // set 3 side pts
         for(ires=1;ires<fit[i].nres-1;ires++){ // New 29jun00
#if 0
           if( (fit[i].nside[ires] == -ARG)||(fit[i].nside[ires] == -LYS) ){
              setside = 2; // set 3 side pts
           }
           else
              setside = 1; //set all side pts
#endif
           setonesidefit(ires,den1,geo,fit[i],geo.set,setside);
         }// ires
//         setsideptsfit(den1,geo,fit[i],geo.set,3); // if 1 set all  if n, set n sidepts
         geo.setval = geo.setmain;
         setptsfit(den1,geo,fit[i],geo.set,0); // set points for main only - not side
         geo.setval = i;
         setendptsfit(den1,geo,fit[i],geo.set); // set points for end of main only - not side
    }//i
     if(iden< firstiden -1) // quit after one lower den then firstiden
//     if(iden< firstiden -2) //New 12dec00  quit after two lower den then firstiden
        goto savebones;
     bone.minden = fminden[iden]; 
     strcpy(bone.bonename,den1.mapname);
     strcat(bone.bonename,".bone"); // bone  name is mapname with .bone_minden appended
     sprintf(ch,"%1.1f",bone.minden);
     strcat(bone.bonename,ch);
     fp=fopen(bone.bonename,"r");
     if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone);
     }
     else{
       fclose(fp);
       readbonefile(bone); // This reopens bone file
     }
//if( (geo.fitn==13)&&(fit[geo.fitn].nres==29) )
//  printf("\n");
     makeeditbones(den1,geo,bone,geo.fitn);
     if(ext.metend>0){
        if(firstiden < 0)
           firstiden = iden;
//if(verbose>=2)
        printf("FOR FIT %d: MET FIT %d  USING DEN %5.3f length = %d\n",
            geo.fitn,ext.metend,bone.minden,ext.length); 
        if(ext.length < 0.7*minlength){ // save if 30% less then previous
           minlength = ext.length;
           bestden = bone.minden;
           *ibestden = iden;
#if 1 // New version - 3july00 - just copie to onemain
           savlength = ext.length;
           savmetend = ext.metend;
           copyoneeditbones(ext.mainone,mainonesav);
#endif
//           copyeditbones(ext,extsav);  // old versin
        }
     } // ext.metend
   } // iden
   savebones:;
   if(minlength < 9000){ // then save one edited bones
#if 1 // New version - 3july00 - just copie to onemain
         ext.length=  savlength;
         ext.metend=  savmetend;
         copyoneeditbones(mainonesav,ext.mainone);
#endif
//      copyeditbones(extsav,ext);
      saveeditbones(maidfile.editbones,testmetend); //this writes only the main chains that are on path to met end      
      printf("FOR FIT %d (nres=%d): MET FIT %d (nres=%d) BEST DEN %5.3f length = %d\n",
            geo.fitn,fit[geo.fitn].nres,ext.metend,fit[ext.metend].nres,bestden,ext.length);
   }
   else{
       printf("   fit %d did not meet fit - exit batchext()\n",geo.fitn);
       *foundedit = 0;  // indicates did not find edited bones
   }
   quitrout:;
}
#if 0
void higherden(struct griddata& den1,struct geometry& geo,int newden)
/*
  See if can fit part of lower den edited bones
   by a higher value.
  1) Starting bone pts  i=2, Use i, i+1 as starting pts in editbones,
      See meet fit. IF so, copy, pts 0,i+1 for original den and
      0, end of new bones to editbones and save.
  2) Do reverse - start at ext.metend, see if can hit fitn.
     If so, reverse bones and save.
*/
{

#endif 

void callbatchextend(struct griddata& den1,
          struct geometry& geo,int *foundedit,int *ibestden)
/*  This is for when batchext is called by itself,
   Not part of extendandfuse()
   geo.fitn is the fit that edit bones are being built for
*/
{
   if(geo.nowset == 0){ // set not already allocated
       imat3(geo.set,den1.griddim[0],den1.griddim[1],den1.griddim[2]);// allocates array 3 dim array (mat)
       geo.nowset = 1;
   }
    readfit(geo,maidfile.fitfile);
    origfitnum = fitnum[0]; 
    if(geo.fitn > fitnum[0]){
        geo.fitn = fitnum[0];
        printf("ERROR: geo.fitn > fitnum[0], setting fitn = fitnum[0]\n");
    }
    batchextend(den1,geo,foundedit,ibestden);
//exit(0);
}

#if 1 
void testconnectionass(int testmet,struct griddata& den1,
       struct geometry& geo,float fusedata[],
       int connum[],int ranki[],int ranki2[],int nterm)
/*
  This version is called only if one of fits is already assined
  nterm = 0 -- cterm extension
  nterm = 1    nterm extension
*/
{ 
   int foundedit,ibestden,maxngap;
   int cnum1;
   int ranki3[MAXRANK+1];
   float rankf[MAXRANK+1];
   int origgapnum;

   origgapnum = geo.gapnum;
   connum[0] = -1;  // default indicates cannot extend
   callbatchextend(den1,geo,&foundedit,&ibestden); //Makes edited bones first call - reads fit,setfit, and extend at geo.fitn
   if(foundedit!=1){
     if(nterm==0){
        if(verbose>=1)
          printf("Assigned fit %d did not meet another cterm fit, setting val2 = -1\n",geo.fitn);
        fit[geo.fitn].val2 = -1;
     }
     else{
        if(verbose>=1)
          printf("Assigned fit %d did not meet another nterm fit, setting val2 = -3\n",geo.fitn);
        fit[geo.fitn].val2 = -3;
     }
     goto quitrout;
   }
   else { //foundedit = 1 found connection
      lowerconnectnumca(geo,ibestden,&cnum1,&maxngap); // should work for both cterm and nterm extension
      if(cnum1<4)
         geo.gapnum = 1;
      else if(cnum1<9)
         geo.gapnum = 2;
      else
         geo.gapnum = 3; 
      if(nterm ==0){ // extend geo.fitn from cterm end to ext.metend
          if(cnum1>maxngap){
             if(verbose>=1)
                printf("SKIPPING assigned fit %d because gap =%d is > %d (setting val2 = -1)\n",geo.fitn,cnum1,maxngap);
             fit[geo.fitn].val2 = -1;
             goto quitrout;
          }
          else{ 
             connectnumass(geo,den1,geo.fitn,ext.metend,connum,ranki,ranki2,ranki3,rankf,
                  ibestden,cnum1,nterm); // if connum (global variable)= gap >0 then has asigned seq to fits
          }
      }
      else{// estend geo.fitn from nterm end back to ext.metend
          if(cnum1>maxngap){
             if(verbose>=1)
                printf("SKIPPING assigned fit %d because gap =%d is > %d (setting val2 = -3)\n",geo.fitn,cnum1,maxngap);
             fit[geo.fitn].val2 = -3;
             goto quitrout;
          }
          else{ 
            connectnumass(geo,den1,ext.metend,geo.fitn,connum,ranki,ranki2,ranki3,rankf,
                  ibestden,cnum1,nterm); // if connum (global variable)= gap >0 then has asigned seq to fits
          }
      }
   } // foundeid
    quitrout:;
   geo.gapnum=origgapnum;
}
#endif


#if 1 
void testconnection(int testmet,struct griddata& den1,
       struct geometry& geo,float fusedata[],
       int connum[],int ranki[],int ranki2[],int nterm)
/*
NOTE - default geo.gapnum is used at all times here
 NEW (30oct00) New allow possibility that fitmet is assigned
    (no longer correct) NEW (25sep00)  THIS IS ONLY CALLED FOR case of metext = 0 (neither fit is assigned)
 New (8jun00) if testmet>0, then if ext.metend does not equal testmet, quit
   Input geo.fitn  = fit # that is being extended in cterm direction
   Returns iconnect = 1 if can definitely assigne seq. # to geo.fitn and
     the fit it meets
   Else returns 0
*/
{ 
   int foundedit,ibestden,maxngap,maxngap2;
   int cnum1,cnum2;
   int fit1,fit2,fit3,ranki3[MAXRANK+1],origusenfits;
   float rankf[MAXRANK+1];

//   maxugap=8; // NEW - now global geo.maxugap   - maximum gap allowed for unassigend cterm extension
   origusenfits = ext.usenfits;
   ext.usenfits=2; // Used for first test
   fit3=0; // NOTE fit3 is not used if ext.usenfits=2;- this is just a dummy value to be used in call
//maxugap =20;
   connum[0] = -1;  // default indicates cannot extend
   callbatchextend(den1,geo,&foundedit,&ibestden); //Makes edited bones first call - reads fit,setfit, and extend at geo.fitn
   if(foundedit!=1){
     if(verbose>=1)
          printf("Unassigned fit %d did not meet another fit, setting val2 = -1\n",geo.fitn);
     fit[geo.fitn].val2 = -1;
     goto quitrout;
   }
   else { //foundedit = 1 found connection
#if 0 //NOT USED - cannot really count on accuracy even if close 9feb01- check if gap connection is 0 or 1 - then use exact gap-
int nres;
float dist;
       if(geo.fixend==2){ // cterm addition      
          nres = fit[ geo.fitn ].nres-1;
          dist = distv( fit[ geo.fitn].pos[nres][1],
                            fit[ ext.metend ].pos[0][1]); 
          if(verbose>=1)
             printf("dist between cterm ends = %5.3f\n",dist);
       }
#endif
       lowerconnectnumca(geo,ibestden,&cnum1,&maxngap); // should work for both cterm and nterm extension
#if 0 // New 30oct00 - try extending even if met assigned fit - 
       if(fit[ext.metend].val1<0){
          if(verbose>=1)
              printf("SKIPPING unassigned fit %d because it met assigned fit %d (setting val2 = -20)\n",geo.fitn,ext.metend);
          fit[geo.fitn].val2 = -20;
          goto quitrout;
       }
       else if(cnum1>geo.maxugap){
#endif
       if(cnum1>geo.maxugap){
          if(verbose>=1)
             printf("SKIPPING unassigned fit %d because gap =%d is >geo.maxugap = %d (setting val2 = -1)\n",
                   geo.fitn,cnum1,geo.maxugap);
          fit[geo.fitn].val2 = -1;
          goto quitrout;
      }
   } 
   // if reach here - test the connection 
   if(fit[ext.metend].val1<0){ // See if unique gap compatible with assignment
      cnum2 = -1; // default, second gap not used
      if(cnum1<4)
         geo.gapnum = 1;
      else if(cnum1<9)
         geo.gapnum = 2;
      else
         geo.gapnum = 3; 
      connectnum(geo,den1,geo.fitn,ext.metend,fit3,connum,ranki,ranki2,ranki3,rankf,
                  ibestden,cnum1,cnum2,nterm); // if connum (global variable)= gap >0 then has asigned seq to fits
      if(connum[0]==1){
         if(verbose>=1)
             printf("Using best connection because fitmet already assigned\n");
      }
      else{ // cannot allow more than 1 connection if fitmet allready assigned 
         if(verbose>=1)
            printf("No unique onnection compatible with fitmet assignment, (setting val2 = -10) for fit%d\n",geo.fitn);
         fit[geo.fitn].val2 = -10;
         connum[0]=0;
      } 
      goto quitrout; // do not try third fit if fitmet assigned
   }//fit[ext.metend].val1<0
   else if( (fit[geo.fitn].nres + fit[ext.metend].nres -4 >=geo.minnres)&&
                  (fit[geo.fitn].nres + fit[ext.metend].nres -4 <=16) ){//met a fit and if at least geo.minnres (13) side chains then try 2 fits 
        cnum2 = -1; // default, second gap not used
        connectnum(geo,den1,geo.fitn,ext.metend,fit3,connum,ranki,ranki2,ranki3,rankf,
                  ibestden,cnum1,cnum2,nterm); // if connum (global variable)= gap >0 then has asigned seq to fits
        if(connum[0]==1){ // found unique fit
           if(ranki2[1]==cnum1 ){ // if  ideal gap  = best gap, then use this result
              if(verbose>=1)
                   printf("Using connection from 2 fits because ideal gap = best gap and nside>=geo.minnres\n");
              goto quitrout;
           }
           else if ((rankf[1]-rankf[2]) >= 1.8*geo.rankgap){ // Very large gap - always use
              if(verbose>=1)
                   printf("Using connection from 2 fits because first much better than second and nside>=geo.minnres\n");
              goto quitrout;
           }
        }
        else if( -rankf[1]>geo.rankmax){ // If fails because too poor best ranking - do not try 3rd fit because this is required conditin to test third
           if(verbose>=1)
             printf("SKIPPING unassigned fit %d best fit has poor ranking (setting val2 = -10)\n",geo.fitn);
          fit[geo.fitn].val2 = -10;
          goto quitrout;
       }
   } // >geo.minnres
   else if( fit[geo.fitn].nres + fit[ext.metend].nres -4 >=17){ // if > 17 - use even if not ideal gap
        cnum2 = -1; // default, second gap not used
        connectnum(geo,den1,geo.fitn,ext.metend,fit3,connum,ranki,ranki2,ranki3,rankf,
                  ibestden,cnum1,cnum2,nterm); // if connum (global variable)= gap >0 then has asigned seq to fits
        if( (connum[0]==1) && ( -rankf[2]> -rankf[1]+ 1.5*geo.rankgap) ){// if > large range between best and second best - allow error in gap if gap>3 (-rank because rank<0)
             if(verbose>=1)
                       printf("Using connection from 2 fits because good rankf and nside>=17\n");
             goto quitrout;
        }
        else if( -rankf[1]>geo.rankmax){ // If fails because too poor best ranking - do not try 3rd fit because this is required conditin to test third
           if(verbose>=1)
             printf("SKIPPING unassigned fit %d best fit has poor ranking (setting val2 = -10)\n",geo.fitn);
           fit[geo.fitn].val2 = -10;
           goto quitrout;
        }
    }
    ext.usenfits=3; // If reach here - need to try 3 fits
    connum[0] = -1;  // Set back to bad default value indicates cannot extend
//exit(0);
    if(verbose>=1)
             printf("Finding third connected fit\n");
    fit1 = geo.fitn;
    fit2 = ext.metend;
    geo.fitn = ext.metend;   
    saveeditbones(maidfile.firstedit, -1); //-1 indicates just save all bone - do not need to find main This will be used later for extension
    callbatchextend(den1,geo,&foundedit,&ibestden); //Makes edited bones first call - reads fit,setfit, and extend at geo.fitn
    if(foundedit!=1){
       if(verbose>=1)
          printf(" Edited bones do not reach third fit, setting val2 = -1 for fit %d  and val2 = -10 for fit %d and quitting\n",fit2,fit1);
       fit[fit2].val2 = -1;
       fit[fit1].val2 = -10;
       goto quitrout;
    }
    else if (ext.metend == fit1){ // comes back on itself
       if(verbose>=1)
          printf(" Edited circle back on itself, setting val2 = -1 for fit %d  and val2 = -10 for fit %d and quitting\n",fit2,fit1);
       fit[fit2].val2 = -1;
       fit[fit1].val2 = -10;
       goto quitrout;
    }
    lowerconnectnumca(geo,ibestden,&cnum2,&maxngap2); // if reach here - found third connection
    fit3 = ext.metend;
    geo.fitn = fit1; // set fitn back to original
    if (cnum2> maxngap2){  // Allow longer gap than 9 here because this gap only used for assignment
             if(verbose>=1)
                  printf(" Edited bones reach third fit but gap = %d > %d: Setting val2=-1 for fit %d- quitting\n",
                                      cnum2,maxngap,fit2);
             fit[fit2].val2 = -1; // Indicates extension failed, do not try again
             savefit();
             goto quitrout;
     }
     else{// third extension successfull
             ext.metend = fit2;
             readeditbones(maidfile.firstedit); // replace bones with original
             saveeditbones(maidfile.editbones, -1); //write original bones to saved file-1 indicates just save all bone - do not need to find main This will be used later for extension
             connectnum(geo,den1,fit1,fit2,fit3,connum,ranki,ranki2,ranki3,rankf,
                  ibestden,cnum1,cnum2,nterm); // if connum (global variable)= gap >0 then has asigned seq to fits   
             if(connum[0]!=1){ // cannot find uniques gap
                if(verbose>=1)
                   printf("Cannot find unique gap for third extension -setting val2 = -10 for fit %d\n",fit1);
                fit[fit1].val2 = -10;
             }
             else if(ranki3[0] == -100){// indicates 3rd fit has sequence assighed, if connum[0]=1, then best gap comparitble with assigment - use this assigment
                 goto quitrout;  //return connum[0]=1
             }
             else if(fit[fit1].nres+fit[fit2].nres+fit[fit3].nres -6 <= geo.minnres){ //(if reach here connum[0]=1)  then ideal gaps must = best gaps
                if( (ranki2[1]!=cnum1)||(ranki3[1]!=cnum2) ){// Only use this fit if ideal gaps = best gaps 
                   connum[0] = -1; // set back to bad value because ideal gaps not = bestgaps
                   if(verbose>=1)
                      printf("Not using this connection because both ideal gaps do not best gaps setting val2 = -10 for fit %d\n",fit1);
                   fit[fit1].val2 = -10;
                }
             }
             else{ // number of sides> geo.minnres   allow 1 gap (>=4) to be non-ideal
                if( (ranki2[1]!=cnum1)&&(ranki3[1]!=cnum2) ){ // both gaps not ideal
                   connum[0] = -1; // set to bad value because both ideal gaps not = bestgaps
                   if(verbose>=1)
                      printf("Not using this connection because neither ideal gaps != best gap, setting val2 = -10 for fit %d\n",fit1);
                   fit[fit1].val2 = -10;
                }
             }                   
    } // third extension successfull
    quitrout:;
    ext.usenfits = origusenfits;
}
#endif

int phiallow(struct geometry& geo,float phi,int restype)
/* Returns 1 if phi in allowed region, 0 if not
 If glycine, any angle allowed
 If proline, -80< phi < -45
 If any other residue - then restricted to alpha helix, beta sheet or
     left helix
*/
{
   int i;
   if(restype == GLY) // any angle allowed
      return 1;
   else if(restype == PRO){
      if( (-80.0*PI/180.0 <= phi)&&(phi<= -45*PI/180.0) )
         return 1;
      else
         return 0;
   }
   else{ // not glycine or proine
     for(i=0;i<=3;i++){
       if( (geo.phiallow[i][0]<=phi) &&(phi<=geo.phiallow[i][1]) )
            return 1;
     }
   }
   return 0; // if reach here, not allowed
}
int psiallow(struct geometry& geo,float phi,float psi,int restype)
/* Returns 1 if psi in allowed region, 0 if not
 If glycine, any angle allowed
 If any other residue - then restricted to alpha helix, beta sheet or
     left helix
*/
{
   int i;
   if(restype == GLY) // any angle allowed
      return 1;
   else{
     for(i=0;i<=3;i++){
       if( (geo.phiallow[i][0]<=phi) &&(phi<=geo.phiallow[i][1]) )
          if( (geo.psiallow[i][0]<=psi)&&(psi<=geo.psiallow[i][1] )  )
             return 1;
     }
   }
   return 0; // if reach here, not allowed
}

void addstrongang(float phi[],float psi[],int itype[],float phist,
        float psist,int *nang)
/* Adds strong angle to nang
*/
{
     float psideg;

     psideg =psist*180.0/PI;
     if(phist<0){
        if( (psideg< 45.0)&&(psideg> -100.0) ){
            *nang = *nang+1;
            phi[*nang] = phist;
            psi[*nang]=  psist;
            itype[*nang]= 1;
        }
        else if(psideg> 45.0 ){
            *nang = *nang+1;
            phi[*nang] = phist;
            psi[*nang]=  psist;
            itype[*nang]= 2;
        }
        else if( (psideg< 111.0)&&(psideg> 9.0 ) ){
            *nang = *nang+1;
            phi[*nang] = phist;
            psi[*nang]=  psist;
            itype[*nang]= 5;
        }
        else if( (psideg< 61.0)&&(psideg> -41.0 ) ){
            *nang = *nang+1;
            phi[*nang] = phist;
            psi[*nang]=  psist;
            itype[*nang]= 6;
        }
     }// phist <0
     else { // phi >0
       if((psideg< 101.0)&&(psideg> -11.0) ){
            *nang = *nang+1;
            phi[*nang] = phist;
            psi[*nang]=  psist;
            itype[*nang]= 3;
        }
        else if(psideg< -91.0 ){
            *nang = *nang+1;
            phi[*nang] = phist;
            psi[*nang]=  psist;
            itype[*nang]= 4;
        }
     } // phi>0
}  

#if 1
void phipsistrong(struct geometry& geo,struct griddata& den1,
      int restype,int *nang,float phi[],float psi[],int itype[])
/* New - 27dec00 - add 2 strongest angles to set of phi/psi previously 
     determine by phipsistart
*/
{
   int i,nstrong,maxstrong;
   float maxphipsi[MAXSTRONG][2], maxden[MAXSTRONG];
   float phideg;

   maxstrong = 2; // maximum number of angles to add to set
   if(geo.fixend==2) // cterm addition
      findstrongphipsicterm(geo,den1,maxphipsi,maxden,&nstrong);
   else 
      findstrongphipsinterm(geo,den1,maxphipsi,maxden,&nstrong);
   for(i=1;i<=nstrong;i++){ // Only angle if maxden big enough
       if(maxden[i]<  0.5*geo.extendden)
          goto badi;
   }
   badi:;
   if(i<maxstrong)
      maxstrong = i;
   for(i=1;i<=maxstrong;i++){
     phideg = maxphipsi[i][0]*180.0/PI;
     if(restype==PRO){   
       if( (phideg> -80.0)&&(phideg< -50.0) )// then add to nang
         addstrongang(phi,psi,itype,maxphipsi[i][0],maxphipsi[i][1],nang);
     }
     else if(restype==GLY){ // then add maxstrong angle
        *nang = *nang+1;
         phi[*nang] = maxphipsi[i][0];
         psi[*nang]=  maxphipsi[i][1];
         itype[*nang]= 0;
     }
     else // Not glycine or Proline             
         addstrongang(phi,psi,itype,maxphipsi[i][0],maxphipsi[i][1],nang);
  }// i
}
#endif 
 
#if 1  // New version
void phipsistart(int restype,int *nang,float phi[],float psi[],int itype[])
/* 
New: use much more limited allowed region for PRO
 Define 2 new types:
   Returns set of allowed phi/psi that can be used
   as input for Powell refinement.
   The points cover allowed regions.
   nang is the number of allowed starting set
   and phi[1..nang] are the angles.
*/
{
   int i,nhelix,nsheet,ndelta,ngly,nleft,neps,ncorhelix,istart;
   float fcv;

   if(restype ==PRO){ //Special case: 1 and 2 helix, 3 and 4 sheet
      phi[1] = -58.0; //Proline helix
      psi[1] = -48.0;
      itype[1]= 8; 
      phi[2] = -58.0;
      psi[2] = -18.0; 
      itype[2]= 8; 
      phi[3] = -58.0;//Proline sheet
      psi[3] = 117.0; 
      itype[3]= 9; 
      phi[4] = -58.0;
      psi[4] = 168.0; 
      itype[4]= 9; 
      *nang =4;
      goto quitrout;
   }
//   nhelix = 2;
 nhelix = 3; //NEW 12oct00 - New add another helix condition 
   ndelta = 2;
   nsheet = 4;
   nleft = 1;
   ngly = 6;
   neps = 2;
   ncorhelix = 1;  // NEW 12oct00 - region at lower right corner of allowed helix region
#if 1 // Used up to 21dec00
   phi[1]= -55.0; // helix
   psi[1]=  -45.0;
   phi[2]=  -110.0;  // helix
   psi[2]=  15.0;
   phi[3]= -130.0;  // NEW 18jan01 - Use different phi 
   psi[3] = -90.0;
#if 0 //Use helix = 4
 nhelix = 4; //NEW 22jan01 - New add another helix condition 
   phi[3]= -60.0;  // NEW 10oct00 - add another helix 
   psi[3] = -90.0;
   phi[4]= -130.0;  // NEW 18jan01 - Use different phi 
   psi[4] = -90.0;
#endif
#endif
#if 0 // tried after 21dec00
   phi[1]= -45.0; // helix
   psi[1]=  -60.0;
   phi[2]=  -110.0;  // helix
   psi[2]=  10.0;
   phi[3]= -130.0;  // NEW 10oct00 - add another helix 
   psi[3] = -70.0;
#endif
    for(i=1;i<=nhelix;i++)
     itype[i]=1;
   istart = nhelix;  
   phi[istart+1]= -70.0;// delta
   psi[istart+1]=  65.0;   
   phi[istart+2]= -120.0;// delta 
   psi[istart+2]=  65.0;   
   for(i=istart+1;i<=istart+ndelta;i++)
     itype[i]=5;
   istart = nhelix+ndelta;  
   phi[istart+1]= -70.0;// sheet
   psi[istart+1]=  115.0;   
   phi[istart+2]= -130.0;// sheet
   psi[istart+2]=  115.0;   
   phi[istart+3]= -70.0;// sheet
   psi[istart+3]=  170.0;   
   phi[istart+4]= -130.0;// sheet
   psi[istart+4]=  170.0;   
   for(i=istart+1;i<=istart+nsheet;i++)
     itype[i]=2;
   istart = nhelix+ndelta+nsheet;  
   phi[istart+1]= 61.0; //left helix
   psi[istart+1]= 45.0;   
   for(i=istart+1;i<=istart+nleft;i++)
     itype[i]=3;
   istart = nhelix+ndelta+nsheet+nleft;  
   phi[istart+1]= 40.0; //epsilon
   psi[istart+1]= -130.0;   
   phi[istart+2]= 90.0; //epsilon
   psi[istart+2]= -179.0;   
   for(i=istart+1;i<=istart+neps;i++)
     itype[i]=4;
   istart = nhelix+ndelta+nsheet+nleft+neps;  
   phi[istart+1]= 0.0; // lower right corner of helix region
   psi[istart+1]= -60.0;  
   for(i=istart+1;i<=istart+ncorhelix;i++)
     itype[i]=10;
   if(restype !=GLY){ 
       *nang = nhelix+nsheet+ndelta+nleft+neps+ncorhelix;
        goto quitrout;
   }
   istart = nhelix+ndelta+nsheet+nleft+neps+ncorhelix;  
   phi[istart+1]= 90.0; //If reach here,must be glycine  The rest of these are for glycine only
   psi[istart+1]= -180.0;   
   phi[istart+2]= 90.0; 
   psi[istart+2]= -135.0;   
   phi[istart+3]= 145.0; 
   psi[istart+3]= -180.0;   
   phi[istart+4]= 120.0; 
   psi[istart+4]= -45.0;   
   phi[istart+5]= 110.0; 
   psi[istart+5]= 160.0;   
   phi[istart+6]= 120.0; 
   psi[istart+6]= 45.0;   
//   for(i=istart+1;i<=istart+ngly;i++) // Used upto 27dec00 
   for(i=1;i<=istart+ngly;i++) // NEW after 27dec00 - set itype = 0 for all angles.
     itype[i]=0;
   *nang = nhelix+nsheet+ndelta+nleft+neps+ncorhelix+ngly; 

   quitrout:;
   fcv = PI/180.0;
   for(i=1;i<= *nang;i++){
      phi[i] = phi[i]*fcv;
      psi[i] = psi[i]*fcv;
   }
}
#endif
   
float caddfunc(float p[],phead)
/*
NO LONGER USED (USED IN POWELL MINIMIZATION)
Funtion to optimize phi, psi for cterm addition
p[0] = delta phi, p[1] = delta psi
*/
{
   float lastden,avecb2,phi,psi,mdist,qual,pos[3],bdist;
   int nbad,nend,mpt[3];

   phi = origphi[0]+ p[1]/5.0;
   psi = origpsi[0]+ p[2]/5.0; 
   makealanine2(geo,geo.pos,0,geo.nend[1],phi,psi);
   geo.nend[1] = geo.nend[1]+1; // Need to set this so avedenreg works OK
   avedenregcb(den1,geo,geo.pos,&lastden,&avecb2,&nbad,
                  geo.nend[1]-1,geo.nend[1]); // Find average MAINCHAIN density for last residue
   geo.nend[1] = geo.nend[1]-1; // Set back to original value
   nend = geo.nend[1];
   maincadist(geo,geo.pos[nend+1][1],geo.main[0],geo.main[1],geo.gnmax,mpt,
                                 &mdist,pos,&bdist); // distance from Ca to main
   if(mdist < 0.7)
      qual= lastden - nbad*maidfile.denscale;
//      qual= lastden;
   else
      qual = -100.0;
//printf("lastden = %5.3f origphi/psi = %5.3f %5.3f, p[1] = %5.3f  p[2] = %5.3f\n",
// lastden,origphi[0],origpsi[0],p[1],p[2]);
   qual =  -qual;//change sign becuse want maximum value, but Powell finds minimum
   return qual;
} 

void comparelastpos(struct geometry& geo,struct pepfit fit1,int fitbeg,
       float *toterr,float *maxerr)
/* Used for cterm extension:
 Compare the main chain pos of last residues in geo and fit
  Returns the total difference in positions
*/
{

   int iend,n0;
   float err;

   iend = geo.nend[1];
   n0 = fitbeg+ geo.nend[1]-geo.nend[0]; // last fit
   *toterr = *maxerr = 0.0;
   err = distv(fit1.pos[n0-1][1],geo.pos[iend-1][1]);// next to last Ca
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0-1][2],geo.pos[iend-1][2]);
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0-1][3],geo.pos[iend-1][3]);
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0][1],geo.pos[iend][1]);// last Ca
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0][4],geo.pos[iend][4]);// last N
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
}

void comparefirstpos(struct geometry& geo,struct pepfit fit1,
       float *toterr,float *maxerr)
/* Used for vterm extension:
  compare the main chain pos of last residues in geo and fit
  Returns the total difference in positions
*/
{

   int iend,n0;
   float err;

   iend = geo.nend[0];
   n0 = 0; // first residue in fit
   *toterr = *maxerr = 0.0;
   err = distv(fit1.pos[n0+1][1],geo.pos[iend+1][1]);// second Ca
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0][2],geo.pos[iend][2]); // C
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0][3],geo.pos[iend][3]); //O
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0][1],geo.pos[iend][1]);// last Ca
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
   err = distv(fit1.pos[n0+1][4],geo.pos[iend+1][4]);// last N
   if(err> *maxerr)
     *maxerr =err;
   *toterr = *toterr + err;
}

void makenewfit(struct geometry& geo,struct pepfit& fitlim,struct pepfit& fits,float qual,float fbadtot,int i)
{
   int j;
          copyfittofit(fitlim,fits,0,fitlim.nres); // copy new fit to fits[i]
          fits.irank[0]=0;  // no previous rankings
          fits.qual=qual;
          fits.fbad= fbadtot;
          fits.fside = 0;
          fits.main[0] = geo.main[0];
          fits.main[1] = geo.main[1];
          fits.main[2] = geo.main[2];
          for(j=1;j<=fitlim.irank[0];j++){ // copy previous ranking and den
             fits.irank[j]=fitlim.irank[j];
             fits.den[j] = fitlim.den[j];
          }
          fits.irank[0]=fitlim.irank[0]+1;// Number of rankings for new addition
          fits.irank[ fits.irank[0] ]= i;  //latest ranking
          fits.den[ fits.irank[0] ] =qual;// latest den
}


void rankfitside(int *nrank,struct pepfit fits[],
        struct geometry& geo,struct griddata& den1,float qualmain,
       float qualside,float fbadtot,
       float mdist,float fitqual[],int endrank,int *totadd)
/*
NEW- 1nov00:  If identical - find aveden, etc for entire fit
    and delete the worse
    totadd = total number of residues added during extension
    ident = 0  - not ident.
    ident = 1  - new is better then identical one
    ident = 2  - new is worse than identical
 ALSO - ASSUME NEW FIT IS INDENTICAL TO AT MOST ONE PREVOUSLY
    RANKED FIT - thus - have commented out prunerank
endrank = total number of previous rankings.
fitqual[1,..,endrank] = previous qualities of fit
New - 1may00:  As go down list of fits, compare main chain atoms pos of last
   residue with each fit.  
   1) Store the positions of this fit.
   2) Go down list in order of highest quality, as before.
   3) If save this fit, continue on down list, comparing positions.
      If find fit with same position at lower point on list, then delete it.
   4) If the fit that currently are comparing quality has same pos, the
         If qual better than fit with same pos, replace it, if worse, then
             delete the new one.
      For this case, do not need to continue down the list
     Note: this may 
Modified from old rankfit()
  - use avedist from bones as criteria for ranking
Try:  weight qualmain and qualside equally.
*/    
{
    int i,ii,j,ident[NAFIT],nres,tadd;
    int prevbad,istart;
    float qual,toterr;
    float limittot,limitmax,maxerr,qual1,qual2;  // max. total and maximum difference in position where consider identical
    float aveden,avecb2,fbad;
    int origmaxside;

    origmaxside = geo.maxside;
    geo.maxside = 3;  // Number of side chain atoms to use when finding density
    limittot =  1.0;
    limitmax = 0.25;
#if 0// For debugging - do not delete similar fits
    limittot =  0.0;
    limitmax = 0.0;
#endif
    tadd = *totadd+1; // number of residues added during extension
    if(geo.fixend==2) // cterm extension
        geo.nend[1] = geo.nend[1]+1;
    else
        geo.nend[0]= geo.nend[0]-1;
#if 0 // Used up to 28sep00
    if(qualside > -1000.0){ // not glycine
//        qual = (qualmain + qualside)/2.0 - 0.5*mdist*maidfile.denscale;// Used up to 3sep00
//        qual = (qualmain + 0.7*qualside)/2.0 - 0.5*mdist*maidfile.denscale; //New decreae weight of qualside
//        qual = (qualmain + qualside)/2.0; // do not mdist
       if(qualside<0) // new 20sep00
           qual = qualmain - 0.5*mdist*maidfile.denscale; //New decreae weight of qualside
       else
           qual = qualmain + 0.7*qualside - 0.5*mdist*maidfile.denscale; //New decreae weight of qualside
//printf("qual = %5.3f\n",qual);    
    }
    else{ 
//        qual = qualmain- 0.5*mdist*maidfile.denscale; 
        qual = qualmain; // do not use mdist
    }
#endif
#if 1 // NEW after 2oct00 - More complicated condition - uses past qual
#if 0
    if(qualside > -1000.0) // not glycine
        qual = qualmain+ 0.7*qualside - 0.5*mdist*maidfile.denscale; //New decreae weight of qualside
    else
        qual = qualmain - 0.5*mdist*maidfile.denscale; //New decreae weight of qualside
#endif
#if 1 // NEW after 31oct00 - Increase weight of mdist from 0.5 to 1.0
    if(qualside > -1000.0) // not glycine
        qual = qualmain+ 0.7*qualside - 1.0*mdist*maidfile.denscale; //New decreae weight of qualside
    else
        qual = qualmain - 1.0*mdist*maidfile.denscale; //New decreae weight of qualside
#endif
    if(qualmain< 0.0){
         if(verbose>=2)
            printf(" Not ranking fit because qualmain =%5.3f < 0.0\n",qualmain);
          goto quitrank;
    }
    else{
        istart = endrank-2; // look at last 3 fits
        if(istart<1)
           istart = 1;
        prevbad =0; // indicates last 3 have good qual
        for(ii=istart;ii<=endrank;ii++){
           if(fitqual[ii]<5.0)
              prevbad =1;// indicates 1 of last 3 have bad qual
        }
        if( (qual<5.0)&&(prevbad==1) ){
           if(verbose>=2)
              printf(" Not ranking fit because qual =%5.3f < 5.0 and previous bad qual\n",qual);
           goto quitrank;
        }
    }
#endif
    *nrank= *nrank+1; // If reach here - rank this fit
    for(i=1;i<= *nrank;i++)
       ident[i] = 0; //default - indicates not identical
    if(geo.fixend==2){ //  cterm extension - copy geo to fitlim
         copyreggeotofit(geo,fitlim,fitbeg); // fitlim (=fit1) was allocated in refineand rank
         fitlim.nres = fitend+2;
    }
    else
         copyreggeotofitnterm(geo,fitlim,fitend); //This also define fitlim.nres
    for(i=1;i< *nrank;i++){ // already increased nrank so no not compare with i = nrank
          if(fits[i].qual <= -1000.0){// indicates has not been allocated to large fit size by copyfit tofit()
             ident[i] = 0;
          }
          else{ 
            if(geo.fixend==2) // cterm extension
               comparelastpos(geo,fits[i],fitbeg,&toterr,&maxerr);
            else
               comparefirstpos(geo,fits[i],&toterr,&maxerr);
//printf("  MAXERR = %5.3f   toterr = %5.3f\n",maxerr,toterr);
             if( (toterr<limittot)&&(maxerr<limitmax) ){ // then identical -only use 1
                ident[i]=1; // this may be changed to 2 below
                goto testident;
             }
           }// else 
     }//i
     testident:;
     for(i=1;i< *nrank;i++){ 
        if(ident[i]==1){    
           if(geo.fixend==2){ // cterm extension
               nres = fitlim.nres-1;
               copyregfittogeo(geo,fitlim,nres-tadd-2,nres);
               avedenregcbside(den1,geo,geo.pos,&aveden,&avecb2,&fbad,
                       geo.nend[1]-tadd,geo.nend[1]); //Returns ave of main and some side density and ave fbad/residue
               qual1 = aveden - 2.0*fbad*maidfile.denscale;
               nres = fits[i].nres-1;
               copyregfittogeo(geo,fits[i],nres-tadd-2,nres);
               avedenregcbside(den1,geo,geo.pos,&aveden,&avecb2,&fbad,
                       geo.nend[1]-tadd,geo.nend[1]); //Returns ave of main and some side density and ave fbad/residue
               qual2 = aveden - 2.0*fbad*maidfile.denscale;
               if(qual1>qual2) // copy fitlim to fits[i]
                   ident[i] = 1; //new fit is better than original
                else
                   ident[i] = 2; // newfit is worse
             } // cterm
             else{ // nterm extension
                copyregfittogeo(geo,fitlim,0,tadd+2);
                avedenregcbside(den1,geo,geo.pos,&aveden,&avecb2,&fbad,
                       geo.nend[0],geo.nend[0]+tadd); //Returns ave of main and some side density and ave fbad/residue
                qual1 = aveden - 2.0*fbad*maidfile.denscale;
                copyregfittogeo(geo,fits[i],0,tadd+2);
                avedenregcbside(den1,geo,geo.pos,&aveden,&avecb2,&fbad,
                       geo.nend[0],geo.nend[0]+tadd); //Returns ave of main and some side density and ave fbad/residue
                qual2 = aveden - 2.0*fbad*maidfile.denscale;
                if(qual1>qual2) // copy fitlim to fits[i]
                      ident[i] = 1; //new fit is better than original
                else
                      ident[i] = 2; // newfit is worse
             } // nterm
             if(ident[i]==2){// newfit is worse -  then just delete this fit
                  *nrank = *nrank-1;
                   goto quitrank;
             }
             else{ // ident = 1 - copy newfit to fit[i]
                  makenewfit(geo,fitlim,fits[i],qual,fbadtot,i);
                  *nrank = *nrank-1;
             }
             goto quitrank;
          }// ident =1
     }//i
     quitident:;
     for(i=1;i<= *nrank;i++){ // if reach here - fits are not identical     
       if( qual>fits[i].qual){ // Insert new fit in this position 
         if(ident[i]==0){// this should always be the case -  bump each previous fit up one if not identical to fit[j]
            for(j= *nrank;j>i;j--){ // move up - this fill up the i=nrank fit
              copyfittofit(fits[j-1],fits[j],0,fits[j-1].nres); // copys fit - including den, qual and nbad data
               fits[j].irank[ fits[j].irank[0] ] = j; // assing last ranking
            }
          } // not identical
          makenewfit(geo,fitlim,fits[i],qual,fbadtot,i);// copy newfit to fits[i]
#if 0
          copyfittofit(fitlim,fits[i],0,fitlim.nres); // copy new fit to fits[i]
          fits[i].irank[0]=0;  // no previous rankings
          fits[i].qual=qual;
          fits[i].fbad= fbadtot;
          fits[i].fside = 0;
          fits[i].main[0] = geo.main[0];
          fits[i].main[1] = geo.main[1];
          fits[i].main[2] = geo.main[2];
          for(j=1;j<=fitlim.irank[0];j++){ // copy previous ranking and den
             fits[i].irank[j]=fitlim.irank[j];
             fits[i].den[j] = fitlim.den[j];
          }
          fits[i].irank[0]=fitlim.irank[0]+1;// Number of rankings for new addition
          fits[i].irank[ fits[i].irank[0] ]= i;  //latest ranking
          fits[i].den[ fits[i].irank[0] ] =qual;// latest den
#endif
#if 0 // Assume only one identity - do not need to prune
          else  // must check if lower ranked fits may be identical
             goto prunerank;
#endif
          goto quitrank; // quit after finding place to insert new fit
       }// if qual > fits[i]
    }// i
#if 0 // Assume only one identity - do not need to prune
    prunerank:;
    for(j=i+2;j<= *nrank;j++){ // start at i+2 because i+1 already compared above
       if(geo.fixend==2) // cterm extension
           comparelastpos(geo,fits[j],fitbeg,&toterr,&maxerr);
       else
           comparefirstpos(geo,fits[j],&toterr,&maxerr);
       if( (toterr<limittot)&&(maxerr<limitmax) ){ // then identical
           ident =1; // indicates identical
           for(k=j+1;k<= *nrank;k++){ // move down
               copyfittofit(fits[k],fits[k-1],0,fits[k].nres); // copys fit - including den, qual and nbad data
           }
           *nrank = *nrank-1;
           goto quitrank;
       }    
    }   
#endif
    quitrank:;
    if(*nrank>= NSAVEFIT)// If nrank over limit, must delete last fit number of saved fits
        *nrank = *nrank-1;
    if(geo.fixend==2) // cterm extension
        geo.nend[1] = geo.nend[1]-1; // set back to value before addition
    else
        geo.nend[0] = geo.nend[0]+1;
    geo.maxside=origmaxside;
}


#if 0
void fixoverlap
/*  Called at beginnig of extend
   Checks all fits for overlap - 
   If overlapped fits meet fuse criteria, then they are fused
       NOTE - at this time - no sequence assignment has been made.
   Else -if cannot fuse, cut end of shorter one until no overlap,
      if shorter one becomes less savlen, then delete it.
   Procedure:  Use routines setendptsfit() and modified version of checkatom
      1) first setptsfit for each fit - setting point to fit number.
      2) then call checkatomlab  - returns 0 if no atom near, and 
           fit number if near fit.
*/
{
   
    for(i=0;i<=fitnum[0];i++){
#endif       



void assignanddyn(struct griddata& den1,struct geometry& geo,int ifit,
     int restype,float *qualmain,float *qualside,float *fbadtot,
     float *mdist,int mpt[],int *lastrank,int *totadd)
/* 
New - 9oct00 - if lastrank <=3, Then use more extensive tordyn at initial step
builds the rotomer and carries of tor. dyn refinement of cterm end
*/
{
   float lastden,avecb2;
   int useset = 1; // means use set when refining
   int nave,nbad,nbadside,resn,bestnbad,bestbadset,bestrotn;
   int totside;
   float pos[3],fbadside;
   float bestqual;
   float aveden,avedenside,firstden;
   float origrfor,bdist,origmaxcbmove,origtorconst,origtorconstside,origtorconstend;
   int origendmove,nstep,ndyn;
   int origmovecb,origverbose;
   int quitcond;
   int icheck;

origverbose= verbose;
//verbose=2;  
//geo.camaindist = 5.0;
   origmovecb = geo.movecb;
   origmaxcbmove = geo.maxcbmove;
   origrfor = tordata.rfor;
   origtorconst = tordata.torconst;
   origtorconstside = tordata.torconstside;
   origtorconstend = tordata.torconstend;
   origendmove = geo.endmove;
   nave = 2; // WORKS WELL: number of resiude to average den over  
//   nave = 3; // number of resiude to average den over
   if(geo.fixend==2){ // cterm extension
      geo.nend[1] = geo.nend[1]+1;
      resn = geo.nend[1]-1; // the side chain that is being built
   }
   else{ // nterm extension
     geo.nend[0] = geo.nend[0]-1;
     resn = geo.nend[0]+1;
   }
   bestrotn = -1; // This means try all rotomers.
   maincadist(geo,geo.pos[geo.nend[1]][1],geo.main[0],geo.main[1],
               geo.gnmax,mpt,mdist,pos,&bdist); // distance from Ca to main
//printf("For geo.rtype = %d: *mdist = %5.3f\n",geo.rtype[resn],*mdist);
#if 0// IT SEEMS IT IS DEFINITELY WORSE TO ALLOW CB TO MOVE!! allow Cb to move in preliminary refineonerotgeo
   geo.movecb = 1; 
   geo.maxcbmove = 0.5;
#endif
   refineonerotgeo(useset,den1,geo,rot,resn,restype,&bestqual,
             &bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
#if 1 // set back to original
   geo.movecb=origmovecb;
   geo.maxcbmove= origmaxcbmove;
#endif
   if(verbose>=2)
//   if(verbose>=1)
            printf("For ifit = %d resn= %d  rottype = %d    bestnbad = %d bestqual = %5.3f\n",
                ifit,resn,restype,bestnbad,bestqual);
//   geo.rtype[resn] =0;// Refine initially with no constraint
   if(*lastrank<=3){  // Use only if not many trial fits
      ndyn = 20;
      geo.endmove = 4; //Definitely better than 2
   }
   else{
      geo.endmove = 2;
      ndyn = 10;
   }
geo.endmove = 4; // New - 24jan01 - all 4 to move initially always - seems to help!!
ndyn = 10;
tordata.rfor= 0.0;
//tordata.torconstside = 0.0; // Turn off side density force during preliminar run
   tordynbat(ndyn,den1,geo,tordata,&aveden,&nbad);//NEW - 4oct00 increase from 10 to 20  Run tordyn- BE CAREFUL quitcond must be >5 if use detalilden below 
   tordata.rfor = origrfor;
   tordata.torconstside = origtorconstside;
   tordata.torconst = origtorconst;
   if(geo.quit==1){ // Something wrong in dynamics - do not use
         if(verbose>=2)// 30nov00 - do not print this lines
             printf("geo.quit =1 in at 1 in assignandyn - setting mdist to bad value, and resetting geo.quit\n");
         *mdist = 1000.0;
         geo.quit = 0; // reset to OK value for next call
    }
    if(geo.fixend==2) // cterm extension
          maincadist(geo,geo.pos[geo.nend[1]][1],geo.main[0],geo.main[1],
               geo.gnmax,mpt,mdist,pos,&bdist); // distance from Ca to main
    else 
          maincadist(geo,geo.pos[geo.nend[0]][1],geo.main[0],geo.main[1],
               geo.gnmax,mpt,mdist,pos,&bdist); // distance from Ca to main
//printf("For geo.rtype = %d: *mdist = %5.3f\n",geo.rtype[resn],*mdist);
   if(*mdist>geo.camaindist){ // minimum condition for using
      *qualmain = -10000.0; // indicates  bad fit, do not rank
       if(verbose>=2)
//       if(verbose>=1)
           printf(" Bad initial mdist = %5.3f\n",*mdist);
   }
   else { // *mdist<=geo.camaindist   -   minimum condition for using
float firstmdist;
       firstmdist = *mdist;
        geo.iside[resn]= -restype; 
        bestrotn = -1; // This means try all rotomers.
        refineonerotgeo(useset,den1,geo,rot,resn,restype,&bestqual,
             &bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
#if 0
//int bestnbad,bestbadset,bestrotn,resn,nbad;
float bestqual,aveden;
 copyreggeotofitnterm(geo,fitlim,fitend);
copyfittofit(fitlim,fit[1],0,fitlim.nres);
savefit();
exit(0);
#endif

#if 1  // try using simulated annealing 
float tormult0;
int i;
//       geo.endmove = 8; 
#if 1// New 29dec00 - do not allow entire trace fit to move - in one case - distortion fromprevious fusion?? cause big shift in positon
         geo.endmove = 4+ *totadd;
         if(geo.endmove>=8)
            geo.endmove=8;
#endif
#if 0 // Used up to July 21, 2000  - works well
       tormult0 = 4.0; // tordata.rfor runs from mult0*rfor to origrfor;
       nstep = 6;
#endif
#if 0 // Used after 28aug00 
       tormult0 = 2.5; //Used 19sep00 - tordata.rfor runs from mult0*rfor to origrfor;
       nstep = 4; // Must be > 1
       quitcond = 6;
#endif
#if 1 // Used after 16jan01 double nstep and halve quitcond (must be >2)
//       tormult0 = 2.5; //Used 19sep00 - tordata.rfor runs from mult0*rfor to origrfor;
// tormult0 = 1.5; //Used 19sep00 - tordata.rfor runs from mult0*rfor to origrfor;
 tormult0 = 2.0; //In one case - definetely better than 1.5
       nstep = 8; // Must be > 1
       quitcond = 3;
#endif
#if 0 //  DEFINITELY WORSE --July 21, 2000  -try higher temp
       tormult0 = 6.0; // tordata.rfor runs from mult0*rfor to origrfor;
       nstep = 8;
#endif
       for(i=1;i<=nstep;i++){
//       tordata.rfor = origrfor*(1.0 + (tormult0-1.0)*(nstep- i)/(nstep-1));
         tordata.rfor = origrfor*tormult0*(nstep- i)/(nstep-1); // runs from 0 to torlult
         tordynbat(quitcond,den1,geo,tordata,&aveden,&nbad);// Run tordyn- BE CAREFUL quitcond must be >2  
         if(geo.quit==1){
             if(verbose>=2)
                printf("geo.quit =1 in at 2 in assignandyn - setting qualmain to bad value, and resetting geo.quit\n");
             *qualmain = -10000.0; // indicates  bad fit, do not rank
             tordata.rfor = origrfor;
             geo.quit = 0; // reset to OK value for next call
             goto endrout;
         }
       }
#endif
      if(geo.fixend==2) // cterm extension
           maincadist(geo,geo.pos[geo.nend[1]][1],geo.main[0],geo.main[1],
                  geo.gnmax,mpt,mdist,pos,&bdist); // distance from Ca to main
       else
           maincadist(geo,geo.pos[geo.nend[0]][1],geo.main[0],geo.main[1],
               geo.gnmax,mpt,mdist,pos,&bdist); // distance from Ca to main
//printf("After sim. anneal: geo.rtype[resn] =%d  mdist = %5.3f\n",geo.rtype[resn],*mdist);
#if 0// SEEMS TO MAKE THING WORSE -  If close to bone, run another round with end force increased
       if(*mdist<geo.camaindist){ // If close to bone, run another round with end force increased
          tordata.rfor = 0;  // Run 1 more round with force increased on end.
          tordata.torconstend = 20*origtorconstend;// New 21dec00 - try increasing density  force on end here
          tordynbat(12,den1,geo,tordata,&aveden,&nbad);// Run tordyn- BE CAREFUL quitcond must be >20 if use detalilden below 
          if(geo.quit==1){
             if(verbose>=1)
                printf("geo.quit =1 in at 2 in assignandyn - setting qualmain to bad value, and resetting geo.quit\n");
             *qualmain = -10000.0; // indicates  bad fit, do not rank
             tordata.rfor = origrfor;
             geo.quit = 0; // reset to OK value for next call
             goto endrout;
          }
        }
        tordata.torconstend=origtorconstend;
#endif
         if(bdist > 7.5){// NEW 4oct00- only check after sim. annealing - prevents from jumping around loop
           if(verbose>=2)
              printf("Bad fit because bdist = %5.3f > 7.5\n",bdist);
           *qualmain = -10000.0; // indicates  bad fit, do not rank
           goto endrout;
       }
       if(*mdist>500.0){// indicates moved to far along bones 
            *qualmain = -10000.0; // indicates  bad fit, do not rank
            goto endrout;
       }
       if(geo.fixend==2) // cterm extension
          avedenregcb(den1,geo,geo.pos,&aveden,&avecb2,&nbad,
                  geo.nend[1]-nave,geo.nend[1]); // Find average MAINCHAIN + Cb density for last tadd residues
       else
          avedenregcb(den1,geo,geo.pos,&aveden,&avecb2,&nbad,
                  geo.nend[0],geo.nend[0]+nave); // Find average MAINCHAIN + Cb density for last tadd residues
       *qualmain= aveden - nbad*maidfile.denscale;
       *qualside = -10000.0;
       *fbadtot = nbad;
       if(restype !=GLY){  //Find side quality for last 2 residues New - fixed call to detaildenoneside on 2may00
          if(*totadd==0){ // then only check side for bad contacts for last residue
             icheck = geo.begend[1]-geo.begend[0]-1; // check for bad contacts for iside>=icheck
             detaildenonesideset(den1,geo,tordata,geo.begend[1]-geo.begend[0]-2,
                geo.begend[1]-geo.begend[0]-1,&avedenside,&nbadside,
                &firstden,&lastden,&totside,icheck); //Correct for both geo.fixend = 1 and 2 finddensity for last(=geo.endmove) side chain - same call used in resqual
//if(nbadside>9)
//  printf("At 1:In assignandyn() nbadside = %d\n",nbadside);
          }
          else{ // check both for bad contacts
             icheck = geo.begend[1]-geo.begend[0]-2;
             detaildenonesideset(den1,geo,tordata,geo.begend[1]-geo.begend[0]-2,
                geo.begend[1]-geo.begend[0]-1,&avedenside,&nbadside,
                &firstden,&lastden,&totside,icheck); //Correct for both geo.fixend = 1 and 2 finddensity for last(=geo.endmove) side chain - same call used in resqual
//if(nbadside>9)
//  printf("AT 2:In assignandyn() nbadside = %d\n",nbadside);
          }
#if 0 // Used up to 15jan01 - the fixend=1 case is wrong - should be identical to fixend=2???
          if(geo.fixend==2) // cterm extension
             detaildenoneside(den1,geo,tordata,geo.begend[1]-geo.begend[0]-2,
                geo.begend[1]-geo.begend[0]-1,&avedenside,&nbadside,
                &firstden,&lastden,&totside); // finddensity for last(=geo.endmove) side chain - same call used in resqual
          else
             detaildenoneside(den1,geo,tordata,1,2,
                &avedenside,&nbadside,
                &firstden,&lastden,&totside); // finddensity for last(=geo.endmove) side chain - same call used in resqual
#endif
          if(totside>=4)
               fbadside = nbadside*4.0/totside;
          else
               fbadside = nbadside;
          *qualside = avedenside - fbadside*maidfile.denscale;
          *fbadtot = nbad + fbadside;
       }
       if(verbose>=2)
//       if(verbose>=1)
          printf("After tordynbat qualmain = %5.3f  qualside = %5.3f fbadtot=%5.3f mdist = %5.3f firstdist = %5.3f bdist = %5.3f\n",
           *qualmain,*qualside,*fbadtot,*mdist,firstmdist,bdist);
#if 0
 copyreggeotofitnterm(geo,fitlim,fitend);
copyfittofit(fitlim,fit[1],0,fitlim.nres);
savefit();
exit(0);
#endif
   }// mdist <1
   endrout:;
   if(geo.fixend==2) // cterm extension
       geo.nend[1] = geo.nend[1]-1;
   else 
       geo.nend[0] = geo.nend[0]+1;
//   geo.maxcbmove=origmaxcbmove;
   quitrout:;
   tordata.rfor=origrfor;
   tordata.torconst=origtorconst;
   tordata.torconstside=origtorconstside;
   tordata.torconstend=origtorconstend;
   geo.endmove=origendmove;
   geo.movecb=origmovecb;
   geo.maxcbmove= origmaxcbmove;
   verbose = origverbose;
}

void addonerescterm(struct griddata& den1,struct geometry& geo,
     struct pepfit fits[],struct pepfit fits2[],int restype,
     int *lastrank,int *nrank,int cis,int metext,int *totadd)
/* routine to add a residue - assume a.a. sequence assigned
New (11apr00) if residue is not a glycine, start with fit # 1 - 
  if nrank >0; then only use fit 1, otherwise, keep trying fits
   until nrank >0
*/
{
   int i,ii,ifit,nend,endrank;
   int nang,mpt[3],main0,main1,resn;
   float phi[MAXSTRONG],psi[MAXSTRONG],mdist;
   int itype[MAXSTRONG];
   float phi0,psi0,fbadtot;
   float qualmain,qualside;
   int k,aanum;  // this line and next line used for cis proline
   float uv[3],vv[3],th,th1,sn,udv;
   float dist1,dist2;
   float fitqual[NAFIT];


//  *nrank = 0;  // NEW  20jun00 - set before call to addoneres
  endrank = fits[1].irank[0];
  for(ifit = 1;ifit<= *lastrank;ifit++){//Try extending each of the saved fits
    for(ii=1;ii<=endrank;ii++)
      fitqual[ii] = fits[ifit].den[ii]; // qual for previous addition
    if(verbose>=2)
      printf("Adding residue at CTERM of fit %d  ( current nrank = %d)\n",
                     ifit,*nrank);
    if(metext!=0) //if = 0 - stop after iconnect added else check enddist
       dist1 = distv(fits[ifit].pos[ fits[ifit].nres-1 ][1],
                   fit[ext.metend].pos[1][1]); //dist for cterm Ca to second Ca in cterm
//printf("dist1 = %5.3f\n",dist1);
    fitend = fits[ifit].nres-1;
    resn = fitend; // the nextside chain residue in fit
    fits[ifit].irot[resn] = 0; // Just rotomer number (only used when save) = 0
    fitbeg = fitend - geo.endmove -3;
    if(fitbeg <0)
       fitbeg = 0;
    phipsistart(restype,&nang,phi,psi,itype);
    copyregfittogeo(geo,fits[ifit],fitbeg,fitend); 
    phipsistrong(geo,den1,restype,&nang,phi,psi,itype); // adds angles that give max den to nang
if(dodebug==1){
   nang=1;
   geo.camaindist = 1000.0;
}
    allocatefit(fitlim,fits[ifit].nres +3); // make large enough so can hold new fit
    main0 = geo.main[0] = fits[ifit].main[0]; // used in maincadist
    main1 = geo.main[1] = fits[ifit].main[1];
//     for(i=1;i<=3;i++){
     for(i=1;i<=nang;i++){
//    for(i=10;i<=nang;i++){
        if(verbose>=2)
//        if(verbose>=1)
           printf("   Using iang = %d  type = %d\n",i,itype[i]);
        copyregfittogeo(geo,fits[ifit],fitbeg,fitend); //NEW
        copyfittofit(fits[ifit],fitlim,0,fits[ifit].nres); // copys fit - including den, qual and nbad data
        nend = geo.nend[1];
        geo.iside[nend]= -restype;
        geo.main[0] = main0; // used in maincadist
        geo.main[1] = main1;
        geo.cis[nend] = cis;
        if(geo.cis[nend]==1){ // if cis proline, rotate last Ca by 180 degrees.
           aanum = nend -1;
           for(k=0;k<=2;k++){
              uv[k]=(geo.pos[aanum][2][k]-geo.pos[aanum][1][k])/geo.dsqrt[1];
              vv[k]=(geo.pos[aanum+1][4][k]-geo.pos[aanum][2][k])/geo.dsqrt[2];
           }
           th= -PI+geo.tt[1];  // standard
//         th1=th+PI-geo.tt[3]; // standard
           th1 = th -PI+ geo.tt[3];  // cis
           udv=cos(th); /* = cn dot nca = cos(theta)*/
           sn=1-udv*udv; /*= sqr(sin(theta)*/
           planevec(uv,vv,udv,sn,geo.dsqrt[4],th,th1,geo.pos[aanum+1][4],geo.pos[aanum+1][1]);
        } // if cis
        makealanine2(geo,geo.pos,0,geo.nend[1],phi[i],psi[i]);
        if(restype!=GLY)
              geo.rtype[nend]=itype[i];
           else
              geo.rtype[nend]=0;
//        powellmin(finalden,pfunc,caddfunc); //(returns best pos = fpos; finalden is pointer Rigid refinemnent returns best tran, rot of ipos
        assignanddyn(den1,geo,ifit,restype,&qualmain,&qualside,
                &fbadtot,&mdist,mpt,lastrank,totadd);
if(dodebug==1){
  qualmain = qualside = 1000.0;
  mdist = 0.0;
}
        if(qualmain > -100.0){  // indicate mdist OK
            if(metext!=0){ //if = 0 - stop after iconnect added else check enddist
              dist2 = distv(geo.pos[geo.nend[1]+1][1],
                   fit[ext.metend].pos[1][1]); //dist for cterm Ca to second Ca in cterm
//printf("dist2 = %5.3f\n",dist2);
              if( (dist1< 4.0)&&(dist2>dist1) ){//This indicate must have extended too far
                if(verbose>=2)
                   printf("Becuuse dist2 (%5.3f) < dist1(%5.3f) - extended too far\n",
                       dist2,dist1);
                goto nexti;
              }
          }
          findphipsi(geo.pos,nend-1,&phi0,&psi0); 
          if(verbose>=2)
              printf("After dyn:phi=%5.3f psi=%5.3f qualmain= %5.3f qualside=%5.3f fbadtot=%5.3f mdist = %5.3f\n",
                  phi0*180.0/PI,psi0*180.0/PI,qualmain,qualside,fbadtot,mdist);

          geo.main[0] = mpt[0]; // these value are stored in fits[] by rankside
          geo.main[1] = mpt[1];
#if 0 // now in rankfitside()
          if(qualmain < 5.0){ // If bad qualmain - discard if previous qual < 5
             if( (endrank>=1)&&(lastqual< 5.0) ){ // NOTE lastqual corresponds to weighed qual, not qualmai
                if(verbose>=1)
                   printf("Bad fit because qualmain = %5.3f <5, and previous qual= %5.3f <5.0\n",qualmain,lastqual);
                goto nexti;
             }
          }
#endif
          rankfitside(nrank,fits2,geo,den1,qualmain,qualside,fbadtot,mdist,
               fitqual,endrank,totadd);
        } // mdist OK
        nexti:;
    }// i
  } // ifit
  savefit:;
}


void addoneresnterm(struct griddata& den1,struct geometry& geo,
     struct pepfit fits[],struct pepfit fits2[],int restype, 
     int *lastrank,int *nrank,int cis,int metext,int *totadd)
/* routine to add residue of type restype at nterm end,
  New (11apr00) if residue is not a glycine, start with fit # 1 - 
  if nrank >0; then only use fit 1, otherwise, keep trying fits
   until nrank >0
*/
{
   int i,ii,ifit,nend,endrank;
   int nang,mpt[3],main0,main1,resn;
   float phi[MAXSTRONG],psi[MAXSTRONG],mdist;
   int itype[MAXSTRONG];
   float phi0,psi0,fbadtot;
   float qualmain,qualside;
   float dist1,dist2;
   float fitqual[NAFIT];

//  *nrank = 0;  // NEW  20jun00 - set before call to addoneres
  endrank = fits[1].irank[0];
//printf("endrank = %d\n",endrank);
  for(ifit = 1;ifit<= *lastrank;ifit++){//Try extending each of the saved fits
    for(ii=1;ii<=endrank;ii++)
      fitqual[ii] = fits[ifit].den[ii]; // qual for previous addition
    if(verbose>=2)
      printf("Adding residue at NTERM of fit %d  (  current nrank = %d)\n",
                     ifit,*nrank);
    if(metext!=0) //if = 0 - stop after iconnect added else check enddist
       dist1 = distv(fits[ifit].pos[0][1],
                   fit[ext.metend].pos[ fit[ext.metend].nres-2 ][1]); //dist for nterm Ca to second Ca in cterm
    fitbeg=0;
    fitend = fitbeg+geo.endmove+3;
    if(fitend> fits[ifit].nres-1)
       fitend = fits[ifit].nres-1;
    resn = fitbeg; // the nextside chain residue in fit
    fits[ifit].irot[resn] = 0; // Just rotomer number (only used when save) = 0
    phipsistart(restype,&nang,phi,psi,itype);
    copyregfittogeo(geo,fits[ifit],fitbeg,fitend); 
    phipsistrong(geo,den1,restype,&nang,phi,psi,itype); // adds angles that give max den to nang
if(dodebug==1){
   nang=1;
   geo.camaindist = 1000.0;
}
    allocatefit(fitlim,fits[ifit].nres +3); // make large enough so can hold new fit
    main0 = geo.main[0] = fits[ifit].main[0]; // used in maincadist
    main1 = geo.main[1] = fits[ifit].main[1];
//for(i=3;i<=3;i++){
    for(i=1;i<=nang;i++){
        if(verbose>=2)
//        if(verbose>=1)
           printf("   Using iang = %d  type = %d\n",i,itype[i]);
        copyregfittogeo(geo,fits[ifit],fitbeg,fitend); //NEW
        copyfittofit(fits[ifit],fitlim,0,fits[ifit].nres); // copys fit - including den, qual and nbad data
        nend = geo.nend[0];
        geo.iside[nend]= -restype;
        geo.main[0] = main0; // used in maincadist
        geo.main[1] = main1;
#if 0
phi[i] = -79.979*PI/180.0;
psi[i] = 118.356*PI/180.0;
#endif
#if 1
        geo.cis[nend] = cis;
//        if(geo.cis[nend]==1) // if cis proline, rotate last Ca by 180 degrees.
#endif
        if(cis==1)
            makecisrevalanine2(geo,geo.pos,0,geo.nend[0],phi[i],psi[i]);
        else
             makerevalanine2(geo,geo.pos,0,geo.nend[0],phi[i],psi[i]);
#if 0
printf("cis = %d\n",cis);
int bestnbad,bestbadset,bestrotn,resn,nbad;
float bestqual,aveden;
geo.nend[0]=geo.nend[0]-1;
resn=geo.nend[0]+1;
refineonerotgeo(1,den1,geo,rot,resn,restype,&bestqual,
             &bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
#if 0
   tordynbat(10,den1,geo,tordata,&aveden,&nbad);// Run tordyn- BE CAREFUL quitcond must be >5 if use detalilden below 
#endif
 copyreggeotofitnterm(geo,fitlim,fitend);
copyfittofit(fitlim,fit[1],0,fitlim.nres);
savefit();
exit(0);
#endif
        if(restype!=GLY)
//              geo.rtype[nend-1]=itype[i];
              geo.rtype[nend]=itype[i];
           else
//              geo.rtype[nend-1]=0;
              geo.rtype[nend]=0;
        assignanddyn(den1,geo,ifit,restype,&qualmain,&qualside,
                &fbadtot,&mdist,mpt,lastrank,totadd);
if(dodebug==1){
  qualmain = qualside = 1000.0;
  mdist = 0.0;
}
        if(qualmain > -100.0){  // indicate mdist OK
          if(metext!=0){ //if = 0 - stop after iconnect added else check enddist
             dist2 = distv(geo.pos[geo.nend[0]-1][1],
                   fit[ext.metend].pos[ fit[ext.metend].nres-2 ][1]); //dist for nterm Ca to second Ca in cterm
             if( (dist1< 4.0)&&(dist2>dist1) ){//This indicate must have extended too far
                if(verbose>=2)
                   printf("Becuuse dist2 (%5.3f) < dist1(%5.3f) - extended too far\n",
                       dist2,dist1);
                goto nexti;
             }
          }
          findphipsi(geo.pos,nend+1,&phi0,&psi0); 
          if(verbose>=2)
//          if(verbose>=1)
              printf("After dyn:phi=%5.3f psi=%5.3f qualmain= %5.3f qualside=%5.3f fbadtot=%5.3f mdist = %5.3f\n",
                  phi0*180.0/PI,psi0*180.0/PI,qualmain,qualside,fbadtot,mdist);

          geo.main[0] = mpt[0]; // these value are stored in fits[] by rankside
          geo.main[1] = mpt[1];
#if 0 // Now in rankfitside
          if(qualmain < 5.0){ // If bad qualmain - discard if previous qual < 5
             if( (endrank>=1)&&(lastqual< 5.0) ){ // NOTE lastqual corresponds to weighed qual, not qualmai
                if(verbose>=1)
                   printf("Bad fit because qualmain = %5.3f <5, and previous qual= %5.3f <5.0\n",qualmain,lastqual);
                goto nexti;
             }
          }
#endif
          rankfitside(nrank,fits2,geo,den1,qualmain,qualside,fbadtot,mdist,fitqual,endrank,totadd);
        } // mdist OK
        nexti:;
    }// i
  } // ifit
  savefit:;
}
void ntermang(float cvec[3],float cavec[3],float posc[3],float posca[3],float *mang)
{
   int k;
   float fitv[3],bonev[3],mag1,mag2;

   for(k=0;k<=2;k++){
     fitv[k] = cavec[k]-cvec[k];
     bonev[k] = posca[k]-posc[k];
   }
   mag1 = sqrt( dot(fitv,fitv) );
   mag2 = sqrt( dot(bonev,bonev) );
   *mang = acos( dot(fitv,bonev)/(mag1*mag2) );
   *mang = *mang*180.0/PI;
}   

void ntermcamainv(float cvec[3],float cavec[3],int pt0,float *ang)
/* Finds angle (in degrees) between nterm C - Ca vec and main bone
  pt0 is is main bone pt at end.
*/
{
   int k;
   float fitv[3],bonev[3],mag1,mag2;

   if(pt0<=0){
     if(verbose>=1)
       printf("WARNING:  mpt = %d <=0 in ntermcamainv()\n",pt0);
     pt0= 1;
   }
   for(k=0;k<=2;k++){
     fitv[k] = cavec[k]-cvec[k];
     bonev[k] = ext.mainone.mainpt[pt0].pos[k] -
                      ext.mainone.mainpt[pt0-1].pos[k];
   }
   mag1 = sqrt( dot(fitv,fitv) );
   mag2 = sqrt( dot(bonev,bonev) );
   *ang = acos( dot(fitv,bonev)/(mag1*mag2) );
   *ang = *ang*180.0/PI;
}   


void addres(struct griddata& den1,struct pepfit fits[NSAVEFIT],
         struct geometry& geo,int *totadd,
         int *lastrank,int metext)
/*
  Add one residue
*/
{
  int i,seq0,resn,restype,cis;
  int nrank;
  float bestqual;
  char rankfile[10],sfit[10],ranknum[10];
  struct pepfit fits2[NSAVEFIT+1];
  struct pepfit fitspro[NSAVEFIT+1];
  int nranktran,mpt;
  float bestqualtran;

   for(i=1;i<NSAVEFIT+1;i++){ //Initial default allocation - automatically reallocated if needed
     fits2[i].nalloc= -10; // initial default - indicates have not yet allocated   
     fits2[i].qual = -1100.0;
     allocatefit(fits2[i],NAFIT);
   }
    seq0 = -fit[geo.fitn].val1;
    if(seq0<=0){
       printf("WARNING:  sequence # not assigned for this fit\n");
       restype = GLY; // if not assigned - allow it to have any phi/psi       
    }
    else{
       if(geo.fixend==2) // cterm addition
          resn = fits[1].nres-1; // all fits should be of same length
       else // nterm
          resn = -(*totadd); // all fits should be of same length
       restype= geo.seq[ seq0 + resn - geo.fnum]; // The added residue typegeo.fnum always  = 1
       if(restype<=0)
          printf("WARNING: a.a. is not assigned for this residue\n");
       else{
         if(verbose>=1)
           printf("Adding residue %d of type %s\n",*totadd+1,nameres[restype]);
       }
    }
    cis = 0; // default use trans
    nrank = 0; // set nrank = 0 before first call to addoneres
    if(geo.fixend==2)
        addonerescterm(den1,geo,fits,fits2,restype,lastrank,&nrank,cis,metext,totadd);// add one residue to all stored fits[]
    else //geo.fixend =1 - nterm
        addoneresnterm(den1,geo,fits,fits2,restype,lastrank,&nrank,cis,metext,totadd);// add one residue to all stored fits[]
#if 0 // Just use for debugging
    if(restype==PRO){ // try cis proline if lastrank <= 2 or bestqual < 20
       printf("PROLINE - nrank = %d\n",nrank);
       for(i=1;i<=nrank;i++){ // Save set of Tran proline fits
          mpt = fits2[i].main[1];
          ntermcamainv(fits2[i].pos[0][2],fits2[i].pos[0][1],mpt,&ang);
          maincadist(geo,fits2[i].pos[0][1],1,mpt-1,geo.gnmax,pt,
                                 &cadist,posca,&ldist); // Ca dist
          maincadist(geo,fits2[i].pos[0][2],1,mpt-3,geo.gnmax,pt,
                                 &cdist,posc,&ldist); // C dist
          ntermang(fits2[i].pos[0][2],fits2[i].pos[0][1],posc,posca,&mang);
          printf("qualilty = %5.3f mpt = %d ang = %5.3f mang = %5.3f cadist = %5.3f cdist = %5.3f\n",
              fits2[i].qual,mpt,ang,mang,cadist,cdist);
       }
    }
#endif
#if 1 // Complicated procedure for trying cis proline
float posca[3],posc[3],ldist,cadist,cdist,cang,maxcdist,pos1[3];
int pt[3],k;
    if(restype==PRO){ // try cis proline if lastrank <= 2 or bestqual < 20
       if(nrank>0){
          bestqual = fits2[1].qual; // the top fit in  fits2.
       }
       else
          bestqual = -1111.0;
       if(verbose>=1)
              printf(" For Proline, best TRAN qual = %5.3f,nrank = %d\n",bestqual,nrank);
       if(geo.fixend==2){ // for cterm extension - use qual as criteria
          if(bestqual<15.0){ // Only  try cis if bad qual
             bestqualtran = bestqual;
             nranktran = nrank;
             if(verbose>=1) // New 16oct00 - always try cis proline.
                 printf(" Trying Cis Prolines\n");
             for(i=1;i<NSAVEFIT+1;i++){ // Allocte fitspro to save set of trans Proline
                fitspro[i].nalloc= -10; // initial default - indicates have not yet allocated   
                fitspro[i].qual = -1100.0;
                allocatefit(fitspro[i],NAFIT);
             }
             for(i=1;i<=nrank;i++){ // Save set of Tran proline fits
                copyfittofit(fits2[i],fitspro[i],0,fits2[i].nres);
                fits2[i].qual = -1100.0; // set fits qual back to 0
             }
             cis=1;
             nrank = 0;  // set back to 0 to count cis prolines
             addonerescterm(den1,geo,fits,fits2,restype,lastrank,&nrank,
                  cis,metext,totadd);// add one residue to all stored fits[]
             if(nrank>0){
                bestqual = fits2[1].qual; // the top fit in  fits2.
             }
             else
                bestqual = -111111.0;
             if(bestqualtran>= 0.8*bestqual){ // use Cis Prolline only if 20% better than tran
               if(verbose>=1)
                   printf(" Using Tran Proline because bestqualtran = %5.3f>=bestqualcis = %5.3f\n",
                         bestqualtran,bestqual);
               bestqual = bestqualtran;
               nrank = nranktran;
               for(i=1;i<=nrank;i++) // Save set of Tran proline fits
                   copyfittofit(fitspro[i],fits2[i],0,fitspro[i].nres);
            } 
            else{ // If cis better, just use current fits2 files     
              if(verbose>=1)
                  printf(" Using Cis Proline because bestqualtran = %5.3f<=bestqualcis = %5.3f\n",
                         bestqualtran,bestqual);
            }
            for(i=1;i<NSAVEFIT+1;i++) 
                freepepfit(fitspro[i]);
         } // bestqual
       } // geo.fixend =2
       else{ //geo.fixend =1 - nterm - more difficult - look at C-Ca angle and mainchain dist and (new) O dist
          if(nrank==0){ // 
            nranktran=nrank;
            bestqualtran = bestqual;
            cang = 1000.0; // set to large value so that will try cis proline
            maxcdist = 1000.0; // These large values will force to use Cis if nrankcis>0 
            goto skipcang;
          }
          mpt = fits2[1].main[1];
          maincadist(geo,fits2[1].pos[0][1],1,mpt-1,geo.gnmax,pt,
                                 &cadist,posca,&ldist); // Ca dist of best fit
          for(k=0;k<=2;k++)
               pos1[k] =ext.mainone.mainpt[ pt[1] ].pos[k]; 
          if( (pos1[0]==posca[0])&&(pos1[1]==posca[1])&&(pos1[2]==posca[2]) ){ // Close bone dist to Ca is not part of bone segment  indicates C-Ca not alligned with bones
              cang = 1000.0; // give it a bad value
              goto skipntermang;
          }
          maincadist(geo,fits2[1].pos[0][2],1,mpt-3,geo.gnmax,pt,
                                 &cdist,posc,&ldist); // C dist
          for(k=0;k<=2;k++)
               pos1[k] =ext.mainone.mainpt[ pt[1] ].pos[k]; 
          if( (pos1[0]==posc[0])&&(pos1[1]==posc[1])&&(pos1[2]==posc[2]) ){ // indicates C-Ca not alligned with bones
              cang = 1000.0; // give it a bad value
              goto skipntermang;
          }
//          maincadist(geo,fits2[1].pos[0][3],1,mpt-3,geo.gnmax,pt,
//                                 &odist,poso,&ldist); // 0 dist
          ntermang(fits2[1].pos[0][2],fits2[1].pos[0][1],posc,posca,&cang);
          skipntermang:;
          if(verbose>=1)
                printf("Tran best qualilty = %5.3f   cang = %5.3f cadist = %5.3f cdist = %5.3f\n",
                  fits2[1].qual,cang,cadist,cdist);
          skipcang:;
//          if( (cang>50.0) || (maxcdist>1.20)||(ocadist< 0.7) ){ // then try cis
#if 0// old version - Use either Cis or Trans, but not both
          if( (cang>50.0) || (maxcdist>1.20)){ // then try cis
             if(verbose>=1) // New 16oct00 - always try cis proline.
                 printf(" Trying Cis Prolines\n");
             bestqualtran = bestqual;
             nranktran = nrank;
             for(i=1;i<NSAVEFIT+1;i++){ // Allocte fitspro to save set of trans Proline
                fitspro[i].nalloc= -10; // initial default - indicates have not yet allocated   
                fitspro[i].qual = -1100.0;
                allocatefit(fitspro[i],NAFIT);
             }
             for(i=1;i<=nrank;i++){ // Save set of Tran proline fits
                copyfittofit(fits2[i],fitspro[i],0,fits2[i].nres);
                fits2[i].qual = -1100.0; // set fits qual back to 0
             }
             cis=1;
             nrank = 0;  // set back to 0 to count cis prolines
             addoneresnterm(den1,geo,fits,fits2,restype,lastrank,&nrank,
                  cis,metext,totadd);// add one residue to all stored fits[]
             if(nrank==0){
                bestqual = -111111.0;
                goto usetran;
             }
             else{ 
               mpt = fits2[1].main[1];
               maincadist(geo,fits2[1].pos[0][1],1,mpt-1,geo.gnmax,pt,
                                 &cadist,posca,&ldist); // Ca dist of best fit
               for(k=0;k<=2;k++)
                 pos1[k] =ext.mainone.mainpt[ pt[1] ].pos[k]; 
               if( (pos1[0]==posca[0])&&(pos1[1]==posca[1])&&(pos1[2]==posca[2]) ){ // indicates C-Ca not alligned with bones
                  nang = 1001.0; // give it a bad value
                  goto skipntermang2;
               }
               maincadist(geo,fits2[1].pos[0][2],1,mpt-3,geo.gnmax,pt,
                                 &cdist,posc,&ldist); // C dist
               for(k=0;k<=2;k++)
                  pos1[k] =ext.mainone.mainpt[ pt[1] ].pos[k]; 
               if( (pos1[0]==posc[0])&&(pos1[1]==posc[1])&&(pos1[2]==posc[2]) ){ // indicates C-Ca not alligned with bones
                 nang = 1001.0; // give it a bad value
                 goto skipntermang2;
               }
               if(cadist>cdist)
                 maxndist = cadist;
               else 
                 maxndist = cdist;
               ntermang(fits2[1].pos[0][2],fits2[1].pos[0][1],posc,posca,&nang);
               skipntermang2:;
               if(verbose>=1)
                  printf("Cis best qualilty = %5.3f  nang = %5.3f cadist = %5.3f cdist = %5.3f\n",
                        fits2[1].qual,nang,cadist,cdist);
               if( (nang<999.0)&&(cang>999.0) )// Use cis if 
                  goto usecis;
               else if( (nang >= cang) || (maxndist >= maxcdist) ){ //Use cis only if both ang and dist are better
                  goto usetran;
               }
               usecis:;
               if(verbose>=1)
                      printf(" Using Cis Proline because better ang and dist\n");
               goto quitpro;
             } // else
             usetran:;
             if(verbose>=1)
                   printf(" Using Tran Proline because better ang or distor nrank=0\n");
             bestqual = bestqualtran;
             nrank = nranktran;
             for(i=1;i<=nrank;i++) // Save set of Tran proline fits
                   copyfittofit(fitspro[i],fits2[i],0,fitspro[i].nres);
             quitpro:; 
             for(i=1;i<NSAVEFIT+1;i++) 
                freepepfit(fitspro[i]);
          } // try cis
#endif
#if 1// New version -18jan01 IF problems with Trans, then use both cis and trans
          if( (cang>50.0) || (maxcdist>1.20)){ // then try cis
             if(verbose>=1) // New 16oct00 - always try cis proline.
                 printf(" Trying Cis Prolines\n");
             cis=1;
             addoneresnterm(den1,geo,fits,fits2,restype,lastrank,&nrank,
                  cis,metext,totadd);// add one residue to all stored fits[]
          } // try cis
#endif
       }// nterm extensio
    }  // Proline
#endif  
    if(nrank>0){ // Then successfully added one  set of nadd residue
#if 0 //Made Worse!! New  - 31oct00 - cut fits that have much poorer quality that the best qual
         bestqual = fits2[1].qual;
         for(i=2;i<=nrank;i++){
            if(fits[i].qual < 0.35*bestqual){
                nrank= i-1;
                goto quitcut;
            }
         }//i
         quitcut:; 
#endif
         geo.extendfitnum = geo.extendfitnum+1; // the set fit # ranked fits are saved in.
         for(i=1;i<=nrank;i++)
            copyfittofit(fits2[i],fits[i],0,fits2[i].nres);
         strcpy(ranknum,geo.setfile);
         strcat(ranknum,"fit");
         strcpy(rankfile,ranknum);
         itoa(geo.extendfitnum,sfit); // sfit is the number of setfile saved in (e.g set2fit(sfit) 
         strcat(rankfile,sfit); 
         saverankfit(rankfile,nrank,fits);
         *lastrank = nrank;  // Number of rankedfits for next iteration
    } //nrank>0 -  savefit
    else{
      *lastrank = 0;
    }
#if 0
printf("QUITTING FOR DEBUGGING PURPOSE\n");
*lastrank= 0;  // this quits routine after one addition
#endif
  for(i=1;i<NSAVEFIT+1;i++) 
      freepepfit(fits2[i]);
    if(*lastrank>0)
       *totadd = *totadd+1;
}


void fusedist(struct pepfit& fit1,struct pepfit& fit2,
       float *maxdist,float *mindist)
/*  
   fit1 is the cterm fit that is being ext, fit2 is the nterm
     fit that was met
returns the maxdist and mindist of atoms that are going to be 
    fused
*/
{
   int i,iend;
   float dist;

   *maxdist =  -1.0;
   *mindist  = 1000.0;
    iend = fit1.nres-1;
    dist = distv(fit1.pos[iend][1],fit2.pos[1][1]); // next to last Ca
    if(dist> *maxdist)
       *maxdist = dist;
    if(dist< *mindist)
       *mindist = dist;
    dist = distv(fit1.pos[iend-1][1],fit2.pos[0][1]); //end Ca
    if(dist> *maxdist)
       *maxdist = dist;
    if(dist< *mindist)
       *mindist = dist;
    dist = distv(fit1.pos[iend][4],fit2.pos[1][4]); // N 
    if(dist> *maxdist)
       *maxdist = dist;
    if(dist< *mindist)
       *mindist = dist;
//    for(i=2;i<=3;i++){ // C, O 
    for(i=2;i<=2;i++){ // C  NEW 10oct00 - do not include O-O dist 
        dist = distv(fit1.pos[iend-1][i],fit2.pos[0][i]); 
        if(dist> *maxdist)
          *maxdist = dist;
        if(dist< *mindist)
          *mindist = dist;
    }
}

void extendassign(struct griddata& den1,struct pepset& fitset,
         struct geometry& geo,int connum0,int metext,int ranki2[],
         int ranki[],float fusedata[],int *totadded,float maxfusedist)
/* 
New (8jun00) returns totadded = #number of residues added
 New routine (15mar00) - extends geo.fitn to ext.metend
     along main bones and makes use of known a.a
  assignment
iconnect = gap number
When have added iconnect+1 residues - fuse the fits (if possible)
*/
{
  int i,lastrank,totadd,bestfit,iconnect,nres;
  float maxdist,mindist,bestdist,bestqual;
  struct pepfit fits[NSAVEFIT];
  int itry; // only 1 fit expected
  int firstadd; // if = 1, then first residue addition
  char fusech[5],fusefile[15];
  float aveden,avecb2,fbad,qual,bestqualall,bestdistall;
  struct pepfit fittemp;
  int iconnectall,totaddall,itry0;
  float tempfusedist,largefusedist,largequal;

//maxfusedist = 4.0;
  largefusedist = 1.7; // fusedist used if good quality
  largequal = -10.0;
  iseed = iseed+1;
  srand48(iseed);// randum number generator (drand48)
  fittemp.nalloc= -10; // initial default - indicates have not yet allocated   
  allocatefit(fittemp,NAFIT);
  if(verbose>=1)
     printf("Starting extendassign using iseed = %d  connum0=%d maxfusedist = %5.3f largefusedist=%5.3f largequal= %5.3f\n",
         iseed,connum0,maxfusedist,largefusedist,largequal);
  for(i=1;i<NSAVEFIT;i++){ //Initial default allocation - automatically reallocated if needed
     fits[i].nalloc= -10; // initial default - indicates have not yet allocated   
     fits[i].qual = -1100.0;
     allocatefit(fits[i],NAFIT);
  }
  copyfittofit(fit[geo.fitn],fits[1],0,fit[geo.fitn].nres); 
  readeditbones(maidfile.editbones); // read the saved edited bones
  fitset.rank[0][0]=0; // number of fits that have associated ranked sets
  ext.ifit = 1;
  fitset.rank[1][0]=0; // No fits in fit 1
  for(i=1;i<NSAVEFIT;i++){ // Set bestqual so  will continue to add new fits
      fits[i].irank[0]=0; // No previous rankings since not a rankedfit
  }
  fits[1].main[0] =1; // only 1 main bone
  fits[1].main[1]=0; // first addition - start on pt 0 of main bone
//  geo.fixend =2; //NOW set before entering fixend = 2 for cterm;= 1 for nterm // refine  cterm end in tordynbat
  totadd = 0;
// following only needed for multcbone
geo.main[0] = fits[1].main[0]; // used in maincadist
geo.main[1] = fits[1].main[1];
  lastrank =1; // This means that for first addition only use fits[1] (the best fit)
  if(metext==0) // then only add iconnect residues and quit
     iconnect = ranki2[1];// If neither fit assigned (metext=0) then add iconnect residues, else - just extend till can fuse - set totadd = iconnect+1
  else // one fit has sequence assinged - extend until can fuse
     iconnect = 10000; // set iconnect to value that can not be satisfied until reset below
  repeat:;
int icnt = 0;
  firstadd = 1;  // indicates that this is first residue added
  bestqualall = -10000.0; // best quality fit that can be fused for all tries
  do{ // If neither fit assigned (metext=0) then add iconnect residues, else - just extend till can fuse - return totadd = iconnect+1
       icnt = icnt+1; 
#if 0 // Use for debugging to look at particular addition in detail
       if(icnt==4)  
           dodebug=1;
#endif
       addres( den1,fits,geo,&totadd,&lastrank,metext);
       bestdist = 1000.0;
       bestqual = -1000.0;
       for(itry=1;itry<=connum0;itry++){
#if 0 // Started writing stuff to cut end if first additin bad, not finished
        if(lastrank>=1) // succesful addtion
            firstadd = 0; // adding second, etc. additon
        else if(firstadd=1){ // first addition failed - cut end and try again
            if(geo.fixend==2){ // cterm addition
               iremove = 2;// remove cterm 
            else 
               iremove = 1;// remove nterm 
            cutend(den1,geo,quittrace,iremove);
        } 
  NEED to add some why of quitting and starting over
#endif
        if(totadd== ranki2[itry]+1){ // check fuse dist
          bestdist = 1000.0;
          for(i=1;i<=lastrank;i++){ // Find fit with least maxdist
            if(geo.fixend==2)// cterm extension
                fusedist(fits[i],fit[ext.metend],&maxdist,&mindist);
            else
               fusedist(fit[ext.metend],fits[i],&maxdist,&mindist);
            if(maxdist < bestdist){
               bestdist = maxdist;
               bestfit = i;
               itry0 = itry;
            }
         }//i
         if(verbose>=1)
             printf("For itry = %d   bestdist = %5.3f\n",itry,bestdist);
        }//totadd = ranki2
      }// itry
      tempfusedist = maxfusedist;
      if(maxfusedist == geo.maxfusedist){// = default = 1.5 - then find bestqual and fuse if good quality and dist < 1.7
          maxfusedist = largefusedist;
      }
      if(bestdist <= maxfusedist){// Then save the bestqual fit
           for(i=1;i<=lastrank;i++){ //New 1nov00  Find fit with best quality
             if(geo.fixend==2){// cterm extension
                fusedist(fits[i],fit[ext.metend],&maxdist,&mindist);
                if(maxdist <=maxfusedist){
                   nres = fits[i].nres-1;
                   copyregfittogeo(geo,fits[i],nres-totadd-2,nres);
                   avedenregcbside(den1,geo,geo.pos,&aveden,&avecb2,&fbad,
                       geo.nend[1]-totadd,geo.nend[1]); //Returns ave of main and some side density and ave fbad/residue
//                   qual = aveden - 2.0*fbad*maidfile.denscale - 1.0*maidfile.denscale*maxdist;
                   qual = aveden - 1.0*fbad*maidfile.denscale - 1.0*maidfile.denscale*maxdist;
                   if(qual>bestqual){
                     bestdist = maxdist;
                     bestfit = i;
                     bestqual = qual;
                   }
                 }
             }// cterm 
             else{ //nterm
                fusedist(fit[ext.metend],fits[i],&maxdist,&mindist);
                if(maxdist <=maxfusedist){
                   copyregfittogeo(geo,fits[i],0,totadd+2);
                   avedenregcbside(den1,geo,geo.pos,&aveden,&avecb2,&fbad,
                       geo.nend[0],geo.nend[0]+totadd); //Returns ave of main and some side density and ave fbad/residue
//                   qual = aveden - 2.0*fbad*maidfile.denscale - 1.0*maidfile.denscale*maxdist;
                   qual = aveden - 1.0*fbad*maidfile.denscale - 1.0*maidfile.denscale*maxdist;
                   if(qual>bestqual){
                     bestdist = maxdist;
                     bestfit = i;
                     bestqual = qual;
                   }
                }        
              }// nterm 
            }// i = lastrank - found fit with best qualith
            if(maxfusedist == largefusedist){ // Case for maxfusedist = 1.7
               if(bestdist>tempfusedist){// tempfusedist = 1.5  - then only fuse if good quality
                  if(bestqual < largequal){
                    if(verbose>=1)
                      printf(" Not fusing because using maxfusedist = %5.3f and bestqual = %5.3f < %5.3f\n",
                               maxfusedist,bestqual,largequal);              
                    goto badqual;
                 }
               } // bestdist >tempfusedist
            } 
            if(bestqual>bestqualall){
               totaddall = totadd;
               bestqualall = bestqual;
               bestdistall = bestdist;
               iconnectall = iconnect;
               copyfittofit(fits[bestfit],fittemp,0,fits[bestfit].nres); // save bestfit
               fittemp.irank[0]=fits[bestfit].irank[0];
               for(i=1;i<=fits[bestfit].irank[0];i++)
                  fittemp.den[i]=fits[bestfit].den[i];
               if(verbose>=1)
                  printf("Extended until fusable,itry=%d lastrank = %d bestdist = %5.3f bestqual= %5.3f\n",
                            itry0,lastrank,bestdist,bestqual);
            }
        } // bestdist<=maxfusedist
        badqual:;
        maxfusedist = tempfusedist;
  } while( (totadd<ranki2[connum0]+1) && (lastrank>0) );
  if(bestqualall> -999.0){ // then have found fusable fit
      *totadded = totadd = totaddall;
      iconnect = iconnectall;
      ranki2[1] = iconnect = totadd-1; // reset ranki2[1] to best gap and this will make jump out of loop.
      if(geo.fixend==2){// cterm extension
              ranki[1] = -fit[geo.fitn].val1+ fit[geo.fitn].nres + iconnect-1;  // Sequence assignment for nterm end of met fit
              if(ranki[1]+fit[ext.metend].nres-2 > geo.maxseq){
                 if(verbose>=1)
                    printf("ERROR: end assignment for fit %d (%d) is >geo.maxseq (%d) -quitting\n",ext.metend,ranki[1]+fit[ext.metend].nres-2,geo.maxseq);
                 lastrank = -1;
              }
       }
       else {// nterm extension
              ranki[1] = -fit[geo.fitn].val1 -fit[ext.metend].nres -iconnect+1;
              if(ranki[1]<geo.seq[0]){
                 if(verbose>=1)
                    printf("ERROR: assignment for fit %d (%d) is <geo.seq[0] (%d) -quitting\n",ext.metend,ranki[1],geo.seq[0]);
                 lastrank = -1;
              }
      }
      if(metext!=0){ // NEW 30oct00 either 1 or both fits previous assigned(if metext=0, only 1 gap allowed and buildandrefine called previouly)
        if(fit[ext.metend].val1<0){ //
           if(verbose>=1)
              printf("Extended until can fuse: Fitmet previously assigned seq = %d \n",ranki[1]);
        }
        else{ // ext.metend not assigned
           if(verbose>=1)
              printf(" Assigning to fitmet:  sequence value = %d   gap = %d\n",
                  ranki[1],ranki2[1]);
           fit[ext.metend].val1 = -ranki[1];
           buildandrefineseq(den1,geo,rot,ext.metend);
           if(geo.quit==1){// out of map range - do not use this extension
              if(verbose>=1){
                  printf(" GOOD EXTENSION BUT NOT USED BECAUSE FIT %d has some side chain out of map range- \n",ext.metend);
                  printf("     reading orig. fit and setting val2 = -3 \n");
              }
              readfit(geo,maidfile.fitfile); //read original fit without this assignment 
              fit[ext.metend].val2 = -3;
              savefit();
              *totadded = totadd;
              fusedata[0] = -1000.0;
              goto quitrout;
           }
        } // ext.metend not assigned
      } // if metext
      countfuse = countfuse+1;
      itoa(countfuse,fusech);
      strcpy(fusefile,"fusefile");
      strcat(fusefile,fusech);
      saverankfit(fusefile,lastrank,fits);
      if(verbose>=1){
         printf("Qual of bestfit = :");
         for(i=1;i<=fittemp.irank[0];i++)
                  printf(" %5.3f",fittemp.den[i]);
         printf("\n");
      }
      printf("Fusing fittemp with %d   maxdist = %5.3f bestqual= %5.3f\n",
                   ext.metend,bestdistall,bestqualall);
      if(geo.fixend==2)// cterm extension
          fuseextfit(den1,geo,0,fittemp,fit[ext.metend],ext.metend);//"0" means remove 0 res from end of bestfit
      else
          fuseextfit(den1,geo,0,fit[ext.metend],fittemp,ext.metend);//"0" means remove 0 res from end of bestfit
      fusedata[0] = bestdist;
  }//if(bestqual>0)
  else{ // bestqualall < -1000 - could not fuse
     if(verbose>=1)
        printf("EXTENSION FAILED BECAUSE bestdist = %5.3f > %5.3f is too large\n",bestdist,maxfusedist);
     *totadded = totadd;
     fusedata[0] = -1000.0;
  }
  quitrout:;
  for(i=1;i<NSAVEFIT;i++) 
      freepepfit(fits[i]);
  freepepfit(fittemp);
}

void savefusedata(struct geometry& geo,int iskip,int fusenum[MAXFUSE][40],
     int *newfuse,int nterm)
/*  newfuse is the new number of the fused fit
*/
{ 
   int i,j,idelete,ifuse,numfitn;
     if( geo.fitn < ext.metend){
            idelete = geo.fitn;
            ifuse = ext.metend;  // the number of saved fused fit
     }
     else{
            idelete = ext.metend;
            ifuse = geo.fitn;
     }
     if(nterm==0)
         numfitn = fusenum[geo.fitn][0];
     else  
         numfitn = fusenum[ext.metend][0];
     if(iskip == 0){ // array not oveflowed
           fusenum[ifuse][0] = fusenum[geo.fitn][0] +fusenum[ext.metend][0]; // number of fits in fused fits
           if(fusenum[ifuse][0] >=39){
              printf("ERROR: fusenum[%d][0]>= 39 - Do not keep track of fusenum anymore\n",ifuse);
              iskip = 1;
              goto skipfusenum;
           }
           if(nterm==0){ //cterm extension
              for(j=fusenum[ifuse][0];j>numfitn;j--) // add metend at c term end
                 fusenum[ifuse][j] = fusenum[ext.metend][j- numfitn ];     
              for(j=1;j<=numfitn;j++)  // add fitn at nterm end
                  fusenum[ifuse][j] = fusenum[geo.fitn][j];
           }
           else { // nterm extension
              for(j=fusenum[ifuse][0];j>numfitn;j--) // add fitn at c term end
                 fusenum[ifuse][j] = fusenum[geo.fitn][j- numfitn ];     
              for(j=1;j<=numfitn;j++)  // add metend at nterm end
                  fusenum[ifuse][j] = fusenum[ext.metend][j];
           }
           for(i=idelete;i<=fitnum[0];i++){ // Copy to lower number value - fitnum[0] has already been reduced in fusefit()
              fusenum[i][0] = fusenum[i+1][0];
              for(j= 1;j<=fusenum[i][0];j++)
                 fusenum[i][j] = fusenum[i+1][j]; 
           }
           *newfuse = ifuse-1;
           printf("Formed new fit # %d:",ifuse-1);
//           fprintf(fpext,"Formed new fit # %d: consisting of original fits:",ifuse-1);
           for(j=1;j<=fusenum[ifuse-1][0];j++){
//              fprintf(fpext,"  %d",fusenum[ifuse-1][j]);
              printf("  %d",fusenum[ifuse-1][j]);
           }
//           fprintf(fpext,"\n");
           printf("\n");
      } // iskip == 0
      skipfusenum:; 
}         

void writesummary(int iskip,int fusenum[MAXFUSE][40])
// Writes list of all current fis
{
   int i,j;
   if(iskip==0){ //Write summary of data if did not overflow arraw - 
//     fpext = fopen(extenddata,"a");    
     for(i=1;i<=fitnum[0];i++){
         printf(" fit # %d (val 1,2 = %d %d)  Corresponds to orignal fits:",i,fit[i].val1,fit[i].val2);
         for(j=1;j<=fusenum[i][0];j++)
            printf("   %d",fusenum[i][j]);
         printf("\n");
     } // i
   } // iskip
}

void extendonefit(int idebug,struct pepset& fitset,struct griddata& den1,
       struct geometry& geo,int fusenum[MAXFUSE][40],float fusedata[],
       int *ifusenum,int iskip)
/* if idebug = 1, quit extending after cterm 
  NEW VERSION - 28sep00:
      1) Redefine val2:  
          1) val2 = -1 indicates either cterm end does not meet another fit
                 or cterm extension failed.
          2) val2 = -2 indicates both cterm and nterm extension failed.
          3) val2 = -10  indicates cterm extension does not meet another fit
                 under conditions when neither fit assigned (use lower
                 gap number limit for this case)
          4) val2 = -20 if unassinged fit meets assigned fit - it is assumed
                that have previously tried and failed to nterm extend other fit 
            (NEW - 30oct00 - now try extending anyway - may succeed in other 
                direction)
     2) Then: each cycle read through all fits -
           1) First look for fits with val1<0 (indicates assigned)
                 if val2>=0 Do cterm extension
                 if val2>= -1  Do nterm extension
                 if val2 = 3, then all done, skip
                    (Actually, this should be always case becasue a given fit
                      should be extended as far as possible in both directions
                      before quitting.)
                    But need other conditions in case start over.
           2) Then look for fits with val1>=0  - not assigned.
              A) if val2>=0: Do batchextend - 
                  1)If it meets a fit that has val1<0 -
                      (already assigned) then quit and set val2 = -20  because this case
                      should have already been tested
                   (NEW - 30oct00 - now try extending anyway - may succeed in other 
                       direction)
                  2) If it meets another fit - but gap too big set val2 = -10
                  3)  else - try to extend - if extension fails, set val2 = -1
              B) If val2 <0 - then quit
           3) This procedure should be completed after one cycle through all fits -
               However - because afer fusion - fits get renumber - Just cycle through
             all fits starting from beginning after each fusion.
  NEW use (definition) of metext:
     metext = 0  - neither fit assigned
     metext =1   - 1 fit assigned.
*/
{
int i,newfuse,maxadded,totadded,isym[NSYM][2];
int origmaxside,metext,testmet,icheck,origfitn;
int ifind,del[5],nsym,iplus;
int connum[MAXRANK+1],ranki[MAXRANK+1],ranki2[MAXRANK+1],nterm;
char ich[5],fusefile[300],str[400];
float brot[NSYM][3][3],btran[NSYM][3];
float maxfusedist;

geo.extendfitnum = 0;  //Now set in extendonefit - This is just the nubmer of set fits are saved in *(if = 1, then e,g. set2fit1 
icheck = 0; // This means do not check files in symfit
maxadded = -1;  // Corresponds to maximum added in previous trys
if(geo.seqnum== -1){ // indicats pdb sequence not yet read - call readsequence
          if(readsequence(geo,maidfile.pdbsequence)!=1){ //read into geo.seq indicates not successfully read
               printf("ERROR: %s cannot be found\n",maidfile.pdbsequence);
               goto quitrout;
          }
}
testmet = -1; // this means do not check the metend - this is not used anymore
ifind = -1;// indicates did fuse at all
do{ // extend geo.fitn at both ends until cannot fuse anymore.
   geo.fnum = 1; // default value = use all the time
   origfitn = geo.fitn;
   connum[0] = -1; // default- do not extend - must set = 1 below if want to extend
   fusedata[0]= -1.0; // value indicates did not fuse.
   if(fit[geo.fitn].val1>=0){ // fit is not assigned
        if(fit[geo.fitn].val2< 0){ // new - val2 = -2 indicates ambigous assignment - go to next fit
           if(verbose>=1)
               printf("SKIPPING unassigned fit %d because val2 = %d\n",geo.fitn,fit[geo.fitn].val2);
           goto quitrout;
        }
        else{ // testconnection return connum[0]=0 if ambigous or connum[0]=1
           geo.fixend=2; // cterm extennsion
           metext = 0; // geo.fitn is not assigned  - tentative indication that neither fit assigned
           nterm = 0; // indicates trying cterm extension
           if(verbose>=1)
               printf("Trying cterm extension of unassigned fit %d\n",geo.fitn);
           testconnection(testmet,den1,geo,fusedata,connum,ranki,ranki2,nterm);// Now only called if fit unassignee
           if(fit[ext.metend].val1<0)  
              metext = 1;
           geo.fixend=2; // cterm extennsion
        } // else              
   } // va1l>0
   else{ // val1<0 - assigned
      metext = 1;   
      if(fit[geo.fitn].val2 == -3){ // Both cterm and nterm failed 
         if(verbose>=1)
            printf("Quitting fit %d because val2 = -3\n",geo.fitn);
      }
      else if(fit[geo.fitn].val2 == -1){ // cterm failed
         nterm = 1;
         geo.fixend =1; // This means extend at N term end
         if(idebug==1)// just quit after trying cterm extension
             goto quitrout;
         if(verbose>=1)
            printf("Trying nterm extension of fit %d\n",geo.fitn);
         testconnectionass(testmet,den1,geo,fusedata,connum,ranki,ranki2,nterm);
      }
      else{// val2 not = -1 or -3 try cterm
         if(verbose>=1)
             printf("Trying cterm extension for fit %d \n",geo.fitn);
         nterm = 0; // indicates trying cterm extension
         geo.fixend=2; // cterm extennsion
         testconnectionass(testmet,den1,geo,fusedata,connum,ranki,ranki2,nterm);
      }            
   }
   if(connum[0]> 0){ //connum[0]  (if>0) here - have connect[0] sequnece to try and extend and fuse
     if(verbose>=1){
         if(metext==0)
           printf("Try to extend using gapnum = %d\n",ranki2[1]);
         else 
           printf("Try to extend usinging tentative  gapnum = %d\n",ranki2[1]);
     }
     if(metext==0 ){ // Neither fit previous assigned - only allowed for case connum[0] = 1
         printf(" Assigned seq. value to fitext = %d   and to fitmet = %d \n",
             ranki[1],ranki[1]+fit[geo.fitn].nres+ranki2[1]-1);
         fit[geo.fitn].val1 = -ranki[1];
         fit[ext.metend].val1 = -(ranki[1]+fit[geo.fitn].nres+ranki2[1]-1);
#if 0
if(  (fit[geo.fitn].val1== -39)&&(fit[ext.metend].val1 == -52) ){
  printf("SETTING VERBOSE = 2\n");
  verbose=2;
}
#endif
         buildandrefineseq(den1,geo,rot,ext.metend);
#if 0
if(verbose==2)
  exit(0);
#endif

         if(geo.quit==1){// out of map range - - reverting back to oldfit - do not use this extension
              if(verbose>=1)
                  printf("  FIT %d has some side chain out of map range -setting val2 = -3\n",ext.metend);
              readfit(geo,maidfile.fitfile); //read original fit without this assignment 
              fit[ext.metend].val2 = -3; 
              savefit();
              goto outofrange;
         }
         buildandrefineseq(den1,geo,rot,geo.fitn);
         if(geo.quit==1){// out of map range - do not use this extension
              if(verbose>=1)
                  printf("  FIT %d has some side chain out of map range -setting val2 = -3\n",geo.fitn);
              readfit(geo,maidfile.fitfile); //read original fit without this assignment 
              fit[geo.fitn].val2 = -3; 
              savefit();
              goto outofrange;
         }
     }
     else if( fit[geo.fitn].val1>=0){ // if metext = 1 and fitext not assigned, assign fit ext
         printf(" Assigned seq. value to fitext = %d \n",
             ranki[1]);
         fit[geo.fitn].val1 = -ranki[1];
         buildandrefineseq(den1,geo,rot,geo.fitn);
         if(geo.quit==1){// out of map range - do not use this extension
              if(verbose>=1)
                  printf("  FIT %d has some side chain out of map range -setting val2 = -3\n",geo.fitn);
              readfit(geo,maidfile.fitfile); //read original fit without this assignment 
              fit[geo.fitn].val2 = -3; 
              savefit();
              goto outofrange;
         }
     }
     origmaxside = geo.maxside;
     geo.maxside=3; //NOTE: MAY WANT TO INCREASE THIS, ESPECIALLY OR SIDE THAT ARE RSIDUES  number of set atoms that are set 
     initsetptsfit(den1,geo);// This sets points for all existing fits
#if 1 // increase allowed fusedist if mgap = 0 or 1
     if(ranki2[1]==0){// if mgap = 0 fuse for any distance
         maxfusedist = 4.0;
     }
     else if(ranki2[1]<=2) // mgap = 1 or 2 
         maxfusedist = 2.0;
#if 0 //19feb01 - keep default fuse dist even if only 1 possible gap
     else if(connum[0]==1) // if only 1 possible gap, allow large dist
        maxfusedist = 2.0;
#endif 
     else // default
         maxfusedist = geo.maxfusedist;
#endif
#if 0
savefit();
exit(0);
#endif
     extendassign(den1,fitset,geo,connum[0],metext,ranki2,ranki,fusedata,
           &totadded,maxfusedist); // extend at cterm or nterm end using assigned sequence info.
     geo.maxside = origmaxside;
     if(totadded>maxadded)
         maxadded = totadded;
     if(fusedata[0]>0){ // fused with another fit - this saves new fits and reads the new fits with new fit num
           ifind = 0;  // indicates at least one fusion - check for sym mates below
           *ifusenum = *ifusenum+1;
           strcpy(str,"copy ");
           strcat(str,maidfile.fitfile);
           strcat(str," ");
           strcpy(fusefile,"fusefit");
           itoa(*ifusenum,ich);
           strcat(fusefile,ich);
           strcat(str,fusefile);
           system(str);  //= "cp maidfile.fitfile fusefit#" 
           if(verbose>=1)
              printf(" Saving current fit file to file %s\n",fusefile);           
           savefusedata(geo,iskip,fusenum,&newfuse,nterm);
           geo.fitn = newfuse;
           writesummary(iskip,fusenum);
           goto fusedfit;
    } // if fused
    readfit(geo,maidfile.fitfile); //if reach here  - extension failed read original fit without this assignment 
    if( nterm==0 ){ // indicate cterm extension failed
        if(verbose>=1) 
             printf("Not using current assignment because cterm extension failed - seting val2 = -1\n");
        fit[geo.fitn].val2 = -1; // if reach here - Indicates tried extension and failed for this fit, do not try again
    }
    else{
        if(verbose>=1) 
             printf("Not using current assignment because nterm extension failed - seting val2 = -3\n");
        fit[geo.fitn].val2 = -3; // if reach here - Indicates tried extension and failed for this fit, do not try again
    }
   } // connum[0]>0
   else{ // connum[0]<=0
     if(metext!=0){ 
        if( nterm==0 ){ // indicate cterm extension failed
          if(verbose>=1) 
             printf("Cannot extend becasue connum[0]=%d<=0 seting val2 = -1\n",connum[0]);
          fit[geo.fitn].val2 = -1; // if reach here - Indicates tried extension and failed for this fit, do not try again
       }
       else{
         if(verbose>=1) 
             printf("Cannot extend becasue connum[0]=%d<=0 seting val2 = -3\n",connum[0]);
         fit[geo.fitn].val2 = -3; // if reach here - Indicates tried extension and failed for this fit, do not try again
       }
     }
   }
   savefit();  // Must save here so val2 is written
   fusedfit:; // jump to here if fuse
   if(idebug==1)// just quit after trying cterm extension
       goto quitrout;
} while( (fusedata[0]>0) || ( (nterm==0)&&(metext!=0) ) ); //  If cterm extension failed and sequence assigned - repeat so try nterm extension
outofrange:;
if(ifind>=0){ // indicates at least one fusion - find symmates
    ifind = 0; // look for all possible sym. mates, do NOT return brot matrix
    iplus = 0; // Means try all ufits, not just one >origfitn
    symfit(iplus,&ifind,icheck,geo,den1,origfitn,isym,&nsym,brot,btran,del);
    for(i=1;i<=isym[0][0];i++){ // found sym mate
          fit[isym[i][0]].val1 = -isym[i][1];
          buildandrefineseq(den1,geo,rot,isym[i][0]);
    }
    if(isym[0][0]>0)
        savefit();
}// ifind
quitrout:; 
}


void extendandfuse(struct pepset& fitset,struct griddata& den1,
       struct geometry& geo,int fusenum[MAXFUSE][40],int *ifusenum,int iskip)
/* 
Main controlling routine:
New: Procedure for finding all possible connections:
  1) Try extending each fit that has sequence assigned.
  2) Then try fits that are not assighned:
      Find connection and the test sequences for gap + or minus 1.
      When test sequence, do not allow values that overlap already
         assigned sequence.
      If can assign sequence - then extend this fit.

First tries each fit - find edited main and test the a.a assigment 
   If, best value is much better than second best - then start extending

  Writes brief summmary to file fusedata
   geo.fitn = number of current fit that is being extended
  When fuse - save fits and readfits with the lower number of 
    of the two fused fits deleted  = idelete
IF succesful fusion leave geo.fitn unchanged
   else increase geo.fitn by 1 andContinue as long 
       as new geo.fitn <= fitnum[0]
 See extendfit.notes (19nov99) for details
fusedata[0] returns information about fusion:
  if >0  then fused and fusedata[0] = fuse error.
  if maxfitdist is too big, set fusedata[0] = -maxfitdist;
  else:
   = -1.0   fprintf(fpext,"  Found main connection to %d (orignum = %d), (length = %d) but extension failed\n",
                                 ext.metend,orignum[geo.fitn],ext.length);
   = -2.0   fprintf(fpext,"   Edited bones do not meet another fit\n");
   = -5.0   fprintf(fpext,"  Found main connection to %d (orignum = %d), but length = %d is too long\n",
                                 ext.metend,orignum[geo.fitn],ext.length);
  Define fusenum[300][40]
    fusenum[i][1,2 ..] = is for the final fit i, the original fits that
     make it up - where fusenum[i][1] is nterm end.
    fusenum[i][0] = number of fits fused in i.
*/
{
   float fusedata[2];
   int i; // extendfit.notes (19nov99) for details
   int initfitnum,initringang,origminiden;
   int idebug,origminres;
   float origrankmax;
   int iexclude,nassign;
   float origrfor,origrankgap;
   char savfile[300];

  origrfor = tordata.rfor;
  initringang = tordata.useringang;
  tordata.useringang = 1;
//tordata.rfor = 1.5*tordata.rfor; // Try varying (Note, side chain of a.a. may be different - see runtordynnew)
printf("STARTING extendandfuse() - Using tordata.rfor = %5.3f * origrfor   tordata.useringang=%d \n",
         tordata.rfor/origrfor,tordata.useringang);
  iexclude =1; // Means only test regions that have not already been assigned when all assignone
  nassign=0; // default value
  idebug = 0; // default=0, if = 1, quit extending after cterm
//idebug = 1; // default=0, if = 1, quit extending after cterm
//idebug=2; // if =2 , set geo.miniden =origiden and do not repeat
#if 0 // for debugging - just run on one fit
//  idebug= 1; // do not try nterm extension
idebug= 0;
  geo.fitn = 24;
  iseed = 259;
  geo.miniden = 0;
  extendonefit(idebug,fitset,den1,geo,fusenum,fusedata,
               ifusenum,iskip);
  goto quitrout;
#endif
   origrankmax = geo.rankmax;
   origminres = geo.minnres; // Mininum number of side chains required to test gap for 2 fits
   origminiden = geo.miniden; // original default value = 0 - corresonds to bone0.8
   geo.miniden = 3; // Minimum value used in looking for connection - 4 corresponds to bone1.4
if(idebug==2){
   printf("WARNING, WARNING WARNING - SETTING idebug=2 in extendandfuse()\n");
   geo.miniden=origminiden;
}
   repeat:;
//goto skipextend; // for debugging
   printf("STARTING EXTEND USING geo.miniden = %d\n",geo.miniden);
      for(i=1;i<=fitnum[0];i++){ // try extending fits with sequence allready assigned
         if( (fit[i].val1<0)&&(fit[i].val2 != -3) ){
            initfitnum = fitnum[0];
            geo.fitn = i;
            extendonefit(idebug,fitset,den1,geo,fusenum,fusedata,
               ifusenum,iskip);
            if(fitnum[0]!=initfitnum) // this indicates have fused fits
               goto repeat;
          }
      } // all assigned fits
      for(i=1;i<=fitnum[0];i++){ // try extending fits with sequences not yet assigned and not problem eexteendign
         if( (fit[i].val1>=0)&&(fit[i].val2 >=0) ){
            initfitnum = fitnum[0];
            geo.fitn = i;
            extendonefit(idebug,fitset,den1,geo,fusenum,fusedata,
               ifusenum,iskip);
            if(fitnum[0]!=initfitnum) // this indicates have fused fits
               goto repeat;
          }
      } // all assigned fits
      writesummary(iskip,fusenum);
      printf("FINISHED - CANNOT EXTEND ANY MORE using geo.rankmax = %5.3f geo.rankgap= %5.3f\n",
                       geo.rankmax,geo.rankgap); // If reach here, finished.
      if(idebug==2){
        goto skiprepeat;
      }
      if( (nassign==0) &&(geo.miniden != origminiden) ){ // If have not previously set miniden = origminiden = 0 ,and not calling after assignone() below - repeat with larger rankmax
         geo.miniden = origminiden; // set to default value = 0 = corresponds to bone0.8
         printf("Resetting all fits: fit.val2 to +1 and repeating with geo.miniden = %d\n",geo.miniden);
         for(i=1;i<=fitnum[0];i++){
                fit[i].val2 = 1;
         }
         savefit();
         goto repeat;
      }
      if( (nassign==0) &&(geo.rankmax <= origrankmax + 0.001) ){ // If have not previous reset rankmax,and not calling after assignone() below - repeat with larger rankmax
//         geo.rankmax = 1.34*geo.rankmax;
         geo.rankmax = 1.5*geo.rankmax;  //11dec00 - try increasing
         geo.minnres = 6; // reduce from original = 13 to 6
         geo.maxugap = 12;  // Increase max. allowed unassighed gap from 8 to 12
         printf("Resetting all fits with fit.val2 = -10 to +1 and repeating with geo.rankmax = %5.3f, geo.maxugap = %d and geo.minnres = %d\n",
              geo.rankmax,geo.maxugap,geo.minnres);
         for(i=1;i<=fitnum[0];i++){
            if(fit[i].val2 == -10)
                fit[i].val2 = 1;
         }
         savefit();
         goto repeat;
      }
      skiprepeat:;
/* AT this point have assigned most fits and all assigned fits have been
    rotated onto unassigned fits which were then assigned - Thus, unassigned
    fits at this point must not overlap the currently assigned fits.
   WARNING - ABOVE IS NOT COMPLETELY CORRECT - IF ROTATED ASSIGNED FIT
     OVERLAPS A FIT BUT DOES NOT FALL IN DENSITY, THEN WILL NOT BE ASSIGNED
*/
skipextend:;
   geo.rankmax = origrankmax;
   geo.minnres = origminres;
   geo.miniden= origminiden;
#if 1 // 12mar01 - No longer all this - better not to assign fits if cannot be connected
   iexclude =1; // Means only test regions that have not already been assigned
   origrankgap = geo.rankgap;
   geo.rankgap = 7.0; // increase for this step from 5.0 to 7.0
   assignone(geo,den1,iexclude,&nassign); //retrurn number of fits assigned
   geo.rankgap = origrankgap;
   if(nassign>0)
      goto repeat;
#endif
   strcpy(savfile,maidfile.fitfile);
   strcat(savfile,"_extend");
   savebackupfit(savfile); // copies fitfile to savfile      
   assignseqment(den1,geo,"extend"); // This tries to assign fits to segments (or subunit)- reorder fits and writes new fit and a .pdb
   quitrout:;
  tordata.rfor=origrfor;
  tordata.useringang=initringang;
}

void extendandfusebatch(struct pepset& fitset,struct griddata& den1,
       struct geometry& geo)
/* 
Main controlling routine:
New: Procedure for finding all possible connections:
  1) Try extending each fit that has sequence assigned.
  2) Then try fits that are not assighned:
      Find connection and the test sequences for gap + or minus 1.
      When test sequence, do not allow values that overlap already
         assigned sequence.
      If can assign sequence - then extend this fit.

First tries each fit - find edited main and test the a.a assigment 
   If, best value is much better than second best - then start extending

  Writes brief summmary to file fusedata
   geo.fitn = number of current fit that is being extended
  extendandfuse returns -1 if not successful, or = metend, the # of fit
     that was met and fused.
  When fuse - save fits and readfits with the lower number of 
    of the two fused fits deleted  = idelete
IF succesful fusion leave geo.fitn unchanged
   else increase geo.fitn by 1 andContinue as long 
       as new geo.fitn <= fitnum[0]
 See extendfit.notes (19nov99) for details
fusedata[0] returns information about fusion:
  if >0  then fused and fusedata[0] = fuse error.
  if maxfitdist is too big, set fusedata[0] = -maxfitdist;
  else:
   = -1.0   fprintf(fpext,"  Found main connection to %d (orignum = %d), (length = %d) but extension failed\n",
                                 ext.metend,orignum[geo.fitn],ext.length);
   = -2.0   fprintf(fpext,"   Edited bones do not meet another fit\n");
   = -5.0   fprintf(fpext,"  Found main connection to %d (orignum = %d), but length = %d is too long\n",
                                 ext.metend,orignum[geo.fitn],ext.length);
  Define fusenum[300][40]
    fusenum[i][1,2 ..] = is for the final fit i, the original fits that
     make it up - where fusenum[i][1] is nterm end.
    fusenum[i][0] = number of fits fused in i.
*/
{
   int i,iskip,fusenum[MAXFUSE][40],ifusenum; // extendfit.notes (19nov99) for details
   int origtorend,origendmove,origmovecb;
   float origpoorrigid;

//geo.rankgap=0; // for debuffing
//geo.rankmax= 100.0;
//iseed = 188; // Used 23feb01 global variable
iseed = 258; // Used upto 27mar01
//iseed= 339;
   fitlim.nalloc = -10;
   ext.usenfits = 2;// Number of unassigned connected that must be linked to assigne seq. value
ext.usenfits = 3;//  I do not think this makes any difference since it reset in other routines -Number of unassigned connected that must be linked to assigne seq. value
   geo.deldist = sqrt( dot(den1.delx,den1.delx) ); // length of den1.delx vector = 1.445 for shad , res = 2.5
   if(geo.deldist < 1.445)
      geo.gnmax = 7*(1.445/geo.deldist) +1;
   else
      geo.gnmax = 7;
   geo.maxugap = 8;  //Maximum gap allowed to connect unassigned fits
   printf("USING geo.maxngap = %d and geo.gnmax = %d  (geo.deldist = %5.3f)\n",geo.maxngap,geo.gnmax,geo.deldist);
   origmovecb = geo.movecb;
   geo.movecb = 0; // do not move Cb in refineonerot.
//     geo.camaindist = 1.5; //WORKS well used upto 9jun00:USE  when use addtillgly -   max dist from main where run dynamics again, etc.
//     geo.camaindist = 2.0; //Used up to 31octooUSE  when use addtillgly -   max dist from main where run dynamics again, etc.
//    geo.camaindist = 1.7; //TRIED USE on 31oct00 - not clear if problem
    geo.camaindist = 1.9; //On 9nov00 - increased from 1.7 to 1.9 - slows things up but probably safer.
//     geo.camaindist = 2.5; //INCREASES TIME, DOES NOT HELP MUCH -New increase allowed distance - 20jul00
//     geo.camaindist = 10.0; //for debugging
srand48(10); // Use same value for each fit for debugging - initialize randum number generator (drand48)
   geo.checkmaindist = 1; //IMPORTANT - in tordynbat - tried to minimie number of Ca and N more than max dist. from edited amin
//   geo.extendfitnum = 0;  //Now set in extendonefit - This is just the nubmer of set fits are saved in *(if = 1, then e,g. set2fit1 
   tordata.useboneconstr = 0;  // if = 1 Use bone constraint on Ca (or N?) in dynamics
//   nadd = 3; //global variable = # of residues added each step 
//   nadd = 2; //global variable = # of residues added each step 
   nadd = 1; //global variable = # of residues added each step 
   origendmove = geo.endmove;
   geo.endmove = 6; // This seems to be better than default vallue (=4)
//   geo.endmove = 3; // Try smaller value
   origtorend = geo.torend;
   geo.torend = 6; // Number of residues whose density is averaged in avectermden
   origpoorrigid = geo.poorrigid;
   geo.poorrigid = 0.25*maidfile.denscale; // Condition for defining geo.npoor
   strcpy(extenddata,maidfile.fitfile);
   strcat(extenddata,"_data");
   fpext = fopen(extenddata,"w");
   fclose(fpext);  // This erases whatever is in extenddata
   readfit(geo,maidfile.fitfile);
#if 0 // 17nov00 - commented this out - is not used since read fits again before saving
   for(i=1;i<=fitnum[0];i++) // set val2 = 0 for all unassinged fits with val2>0
     if(fit[i].val1>=0)
        if(fit[i].val2>0)
           fit[i].val2=0;  
#endif
#if 0 // for degugging  of fusefit
   origfitnum = fitnum[0];
   geo.fitn = 5;
   fuseextfit(den1,geo,0,fit[5],fit[9],9);
   goto quitrout; 
#endif
   if(fitnum[0]> MAXFUSE-2){ //This conditiion is just needed for orignum[]
      printf("ERROR: fitnum[0]>MAXFUSE-2, CANNOT KEEP TRACK DATA ABOUT WHICH FITS ARE FUSED\n");
      iskip = 1;
   }
   else{
     iskip = 0;
     for(i=1;i<=fitnum[0];i++){
        fusenum[i][0] = 1; // just one fit in each fit (no fusions yet)
        fusenum[i][1] = i;
     }
  }
   geo.fixend =2; // This means extend at c term end
   ifusenum = 0; //number of saved fusefit file
   extendandfuse(fitset,den1,geo,fusenum,&ifusenum,iskip);
   geo.poorrigid = origpoorrigid;
   geo.torend = origtorend;
   geo.endmove=origendmove;
   geo.movecb = origmovecb;
   quitrout:;
}
  
