/* extende.c++    Extra routines used for adding several residues
    at a time.
  Routines that were in old extendfit.c++ (extendfit.c25jan00)
*/


#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 "initialize.h"
#include "bone.h"
#include "rankside.h"
#include "assignseq.h"
#include <string.h>
#define NADDTYPE 100  // max number of saved - classified sets of 3 additions
#define MAXADD 4 // max. number residues added each step
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 int origfitnum;  //keeps track of of numer of fit in fit file
extern int nadd; // global variable  = # of residues to to add at each step
extern int fitbeg, fitend; // = beg and end values in copyregfittogeo,e tc.
extern int totadd;  // total # of res. added after reading original fit 
extern int startres;// Global variable - fit res # where start extenting
extern char extenddata[200];
extern struct pepfit fitlim; //global - becuase use in maxfitdistgeo() called from tordyn
extern int prevbestiden;
extern int iextend; // for degugging

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];
};

#if 1
void checkdist(struct geometry& geo,int nmax,int testpt[],float *avedist)
/*  Check the N and Ca dist of nadd residues afterrefinment
  Returns 0 if not closeenough
 ASSUMES ONLY 1 MAIN EDITED BONE = 1
New - use tespt[0] for first residue
  and new value mpt[ii][1] for subsequent residues
*/
{
    int ii,rend,mpt[MAXADD][3];
    float cadist,ndist,pos[3],bdist;
//    float limitdist = 1.0;
//    limitdist = geo.camaindist;

    *avedist = 0.0;
    rend = geo.nend[1]-nadd+1;
    for(ii = 0;ii<nadd;ii++){ // check last added residues 
       maincadist(geo,geo.pos[rend+ii][1],1,testpt[ii],nmax,mpt[ii],
                                 &cadist,pos,&bdist); // Ca dist
       *avedist = *avedist +cadist;
//printf("Ca dist = %5.3f\n",cadist);
#if 0
       if(cadist>limitdist)
          return 0;
#endif
       maincadist(geo,geo.pos[rend+ii][4],1,testpt[ii],nmax,mpt[ii],
                                 &ndist,pos,&bdist); // N dist
//printf(" N dist = %5.3f\n",ndist);
       *avedist = *avedist +cadist;
#if 0
       if(ndist>limitdist)
          return 0;
#endif
       testpt[ii+1] = mpt[ii][1]; //  Value to use for next residue
    }
    *avedist = *avedist/(2.0*nadd);
}
#endif

void rankfit(int *nrank,struct pepfit fit1,struct pepfit fits[],
        float lastden,int nbad,struct geometry& geo,int val[2],
        float avedist)
/*
New - use avedist from bones as criteria for ranking
New - rank by qual rather then aveden
New - include the main bone (imain) that fit is on
New - allow at most 3 saved fits have cbden = 0 
This allows for more variey in saved fits
  Ranks the set of fits:
  Input latest geo, and its qual, and then
    goes through list of nrank saved fits,and their corresponds quality.
  bestqual[1] is bestfit, bestqual[2] is next best, etc.
nrank =  the number of prevous saved fits, upto a maximum of (e.g. NSAVEFIT-1
qualside includes ave length of side chain in ranking
   seems to make things worse.
*/    
{
    int i,j;
    float qual;

    qual = lastden*(1.0 - ( (float)2.0*nbad/nadd ) - avedist/2.0);
//    qual = lastden*(1.0 - ( (float)2.0*nbad/nadd ));
    for(i=1;i<= *nrank;i++){
       if( qual>fits[i].qual){
//          copyfittofit(fit1,fittemp,0,fit1.nres);
          copyreggeotofit(geo,fitlim,fitbeg); // fitlim (=fit1) was allocated in refineand rank
          fitlim.nres = fitend+nadd+1;
          for(j= *nrank;j>i;j--){
             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
          }
          copyfittofit(fitlim,fits[i],0,fitlim.nres);
          fits[i].val1 = fits[i].val2 = 0;
          fits[i].irank[0]=0;  // no previous rankings
          fits[i].qual=qual;
          fits[i].fbad= nbad;
          fits[i].fside = 0;
          fits[i].cbden = geo.cbden[ geo.nend[1]-1 ];
          fits[i].main[0] = geo.main[0];
          fits[i].main[1] = geo.main[1];
          fits[i].main[2] = geo.main[2];
#if 0
int n0= geo.nend[1];
printf("n0=%d  geo.rtype[n0-2]=%d geo.rtype[n0-1]=%d geo.rtype[n0]=%d\n",
   n0, geo.rtype[n0-2],geo.rtype[n0-1],geo.rtype[n0]);
#endif
//NEW          copygeotofit(geo,fits[i],val); // Save best fit (also set irank[0]=0)
          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] ] =lastden;// latest den
          goto quitrank;
       }// if qual
    }// i
    quitrank:;
}

#if 1 // standard version
int classifyang(struct geometry& geo,float phi,float psi)         
/*
New - try not using type 5:
  and possibly returning more 0 region
 classify phi/psi into constaint types
  Assume -PI <= ang <= PI
  Use large crude range range - may be more than one type
  if(phi<0, and not type 1 or 5, then = type 2 (sheet)
  Types 1,2 and 5 nearly cover phi<0 range
*/
{
#if 0 // use type 5
   if(phi<=0){
     if( (psi<= 40.0*PI/180) &&(psi>= -90*PI/180) )
        return 1;
     else if( (psi> 40.0*PI/180) &&(psi< 85*PI/180) )
        return 5;
     else
        return 2;
   }
#endif
#if 1 // do not use type 5
   if(phi<=0){
     if( (psi<= 65.0*PI/180) &&(psi>= -90*PI/180) )
//     if( (psi<= 45.0*PI/180) &&(psi>= -100*PI/180) )
        return 1;
     else
        return 2;
   }
#endif
   else if( (phi> 20.0*PI/180)&&(phi<100.0*PI/180) ){
      if( (psi<100.0*PI/180) &&(psi> -40*PI/180) )
         return 3;
      else if (psi< -100.0*PI/180)
         return 4;
   }
   return 0; // not any above constraint
}  
#endif


void sortbest(struct pepfit fit1,float fpos[NAFIT][NAA][3],
    struct griddata& den1,struct geometry& geo,
    struct besttype best[NADDTYPE],float phi[MAXADD],float psi[MAXADD],
    int mpt[MAXADD][3],int *iquit,float avedist)
/* 
New: avedist = average dist from edited bone - use as sort parameter
for each geo - store fit with best density for
    each constraint thpe
*/
{
  int i1,i,n,nbad,type[MAXADD];
  float aveden,avecb2,qual;

   *iquit = 0; // if = 1, filled up set, quit adding
   for(i=0;i<nadd;i++){
       if(phi[i]>PI)
          phi[i] = phi[i]-2.0*PI;
       if(psi[i]>PI)
          psi[i] = psi[i]-2.0*PI;
   }//i
   int n0 = geo.nend[1]+nadd-1; 
   for(i=0;i<=2;i++){ // copy last psi rotation to geo
     geo.pos[n0][3][i]=fpos[n0][3][i];
     geo.pos[n0+1][4][i]=fpos[n0+1][4][i];
     geo.pos[n0+1][1][i]=fpos[n0+1][1][i];
   }
#if 0 // This is identical to result when commented out, about 5% slower with it
   copyfittogeo(geo,fit1);
   for(i=0;i<nadd;i++) // each addition
      makealanine2(geo,geo.pos,0,geo.nend[1]+i,phi[i],psi[i]);
#endif
#if 1 // if = 1 use aveden and nbad as sort paramter 
   geo.nend[1]=geo.nend[1]+nadd;
   avedenreg(den1,geo,geo.pos,&aveden,&avecb2,&nbad,
           geo.nend[1]-nadd,geo.nend[1]); // Find average density for last nadd additions
//   qual = aveden  - ( (float)2.0*nbad/nadd )*aveden;
   qual = aveden;
   geo.nend[1] = geo.nend[1]-nadd;
   if(qual>0){// minimum required for refinement
#endif
#if 0 // use avedist as sort paramter - DOES NOT HELP
   if(avedist<= geo.camaindist - 0.2){
       qual = -avedist;
#endif
       for(i=0;i<nadd;i++){ // each addition
           type[i]=classifyang(geo,phi[i],psi[i]);
//           printf("for iadd = %d  type = %d phi = %5.3f  psi = %5.3f\n",i,type[i],phi[i]*180/PI,psi[i]*180/PI); 
       }
       for(i1=1;i1<=best[0].main[0];i1++){//each stored type        
             for(i=0;i<nadd;i++){ // each addition
                if(type[i]!=best[i1].type[i])// New type, add if different from rest of best
                   goto nexti1;
             }//i
             if(qual>best[i1].qual){ // if reach here, same as previous stored type
             	best[i1].qual = qual;
             	best[i1].nbad = nbad;
             	for(i=0;i<nadd;i++){
                   best[i1].phi[i]= phi[i];
                   best[i1].psi[i]= psi[i];
                }//i
                for(i=0;i<3;i++)
                  best[i1].main[i]=mpt[nadd-1][i];
            }//qual>best[i1].qual
            goto nextiset;
            nexti1:;
       }//i1
       n = best[0].main[0] = best[0].main[0] +1; // If reach, here, this is new type
       best[n].qual = qual;
       best[n].nbad = nbad;
       for(i=0;i<nadd;i++){
               best[n].phi[i]= phi[i];
               best[n].psi[i]= psi[i];
               best[n].type[i]=type[i];
       }// i
       for(i=0;i<3;i++)
              best[n].main[i]=mpt[nadd-1][i];
       if(n>= NADDTYPE-1){
               printf("ERROR: IN sortbest, too many types n= %d\n",n);
               *iquit=1;
               goto quitrout;
        }
        nextiset:;
    }//qual
    quitrout:;
}
void rankfitside2(int *nrank,struct pepfit fits[],
        struct geometry& geo,float qualmain,float qualside,float fbadtot,
       float mdist)
/*
Slightly modified versin of rankfitside to be used with multbone
*/    
{
    int i,j;
    float qual;

    if(qualside > -1000.0){ // not glycine
        qual = (qualmain + qualside)/2.0 - 0.5*mdist*maidfile.denscale;
//        qual = (qualmain + qualside)/2.0; // do not mdist
    }
    else{ 
//        qual = qualmain- 0.5*mdist*maidfile.denscale; 
        qual = qualmain; // do not use mdist
    }
    if(qual <5.0){ // do not rank
        if(verbose>=1)
           printf(" Not ranking fit because qual =%5.3f < 5\n",qual);
        *nrank = *nrank-1; // remove last rank added before call
        goto quitrank;
    }
    for(i=1;i<= *nrank;i++){
       if( qual>fits[i].qual){
          copyreggeotofit(geo,fitlim,fitbeg); // fitlim (=fit1) was allocated in refineand rank
          fitlim.nres = fitend+1+nadd;
          for(j= *nrank;j>i;j--){
             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
          }
          copyfittofit(fitlim,fits[i],0,fitlim.nres);
//          fits[i].val1 = fits[i].val2 = 0;
          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
          goto quitrank;
       }// if qual
    }// i
    quitrank:;
}

#if 1 // Version that tests Ca, and n dist 
void addnresidue(struct griddata& den1,struct geometry& geo,struct pepfit fit1,
     int imain,int pt0,int nmax,int nrot,
     struct besttype best[NADDTYPE],int testpt[])
/*
   Add  n=3 residues to geo, and try all nrot phi/psi angles, if endca is
       withing geo.maincadist of mainbone, then add another residue, etc.
   Returns:  nset[0][0][0] = number in set
             nset[i=1,...nset][0,1,2][0,1,2] = main number and point number and (possible) side num
                  of closest mainbone pt for 0 ith residue, 1 residue, 2 rsidue
             phipsi[i=1,...nset][0,1,2][0,1] = phi (0) and psi (1)
*/
{
   int i,ii,nh[MAXADD],np[MAXADD],k,mpt[MAXADD][3],nend[MAXADD];
   float delrot[MAXADD],phi0[MAXADD],psi0[MAXADD],phi[MAXADD],psi[MAXADD],rotv[MAXADD][3],mag;
   float pos[3],mdist,ccdist,nndist;
   float fpos[NAFIT][NAA][3];
   int mrot[MAXADD],iquit;
   float origmaindist,bdist;
   int tmain,tpt;

origmaindist = geo.camaindist;
repeat:;  // goto here if repeat with larger geo.camaindist 
clock_t time0,time1;
   best[0].main[0]=0; // this = total number added (1,2,..n)
   for(i=0;i<NADDTYPE;i++)
       best[i].qual = -1000.0;
   nend[0]=geo.nend[1];
   nend[1]=geo.nend[1]+1;
   nend[2]=geo.nend[1]+2;
   nend[3]=geo.nend[1]+3;
#if 1  // Works well in most cases
   mrot[0]= 2*nrot;
   mrot[1]= 1.5*nrot;
   mrot[2]= nrot;
   mrot[3]= nrot;
#endif
#if 0  //Increase number of rotations - Does not seem to help
   mrot[0]= 3*nrot;
   mrot[1]= 2.0*nrot;
   mrot[2]= 1.5*nrot;
   mrot[3]= nrot;
#endif
  if(nadd <1)
      goto quitrout;
   if(nadd>3)
      printf("ERROR: nadd = %d >3, setting = 3\n",nadd);
   for(i=0;i<nadd;i++)
     delrot[i] = 2.0*PI/mrot[i];
   ii=0; // first step 
   phi0[ii]=psi0[ii] = 0.0; // set phi psi = 0
time0 = clock();

float mtempmin[5],ntempmin[5];
for(i=1;i<=4;i++)
   mtempmin[i]=ntempmin[i] = 1000.0;

   for(nh[ii] =0;nh[ii]<mrot[ii];nh[ii]++){// rotate phi in steps of 10 degrees - find max C den and corresponding Cb den
if (verbose>=2)
     printf("nh[%d] = %d\n",ii,nh[ii]);

     phi[ii]= phi0[ii]+delrot[ii]*nh[ii]; 
     makealanine2(geo,geo.pos,0,nend[ii],phi[ii],psi0[ii]);
     for(k=0;k<=2;k++) // Ca to C vector
        rotv[ii][k] = geo.pos[ nend[ii] ][2][k]-geo.pos[ nend[ii]][1][k];
     mag = sqrt( dot(rotv[ii],rotv[ii]) );  
     for(k=0;k<=2;k++) 
          rotv[ii][k]=rotv[ii][k]/mag;
     for(np[ii]=0;np[ii]<mrot[ii];np[ii]++){
        psi[ii]= psi0[ii]+delrot[ii]*np[ii];
        rotonca(geo,geo.pos,fpos,nend[ii],rotv[ii],geo.pos[nend[ii]][2],psi[ii]); //rotate o,n,and ca about rotv
        testpt[0] = pt0;
        maincadist(geo,fpos[nend[ii]+1][4],imain,pt0,nmax,mpt[ii],
                                 &nndist,pos,&bdist);
ccdist = 0.0;
        maincadist(geo,fpos[nend[ii]+1][1],imain,pt0,nmax,mpt[ii],
                                 &mdist,pos,&bdist);
if(mdist< mtempmin[1]){
   mtempmin[1]= mdist;
}
if(nndist< ntempmin[1]){
   ntempmin[1]= nndist;
}
        if( (mdist< geo.camaindist)&&(ccdist< geo.camaindist)&&
              (nndist< geo.camaindist)  )   { //  test another residue
          if(nadd==1){
             sortbest(fit1,fpos,den1,geo,best,phi,psi,mpt,&iquit,
                   (mdist+nndist)/2.0);
             if(iquit==1)
               goto quitrout;
             goto nextii1;
          }
           makealanine2(geo,geo.pos,0,nend[ii],phi[ii],psi[ii]); //if reach here add good residue
           ii=1; // first step 
           phi0[ii]=psi0[ii] = 0.0; // set phi psi = 0
           for(nh[ii] =0;nh[ii]<mrot[ii];nh[ii]++){// rotate phi in steps of 10 degrees - find max C den and corresponding Cb den
             phi[ii]= phi0[ii]+delrot[ii]*nh[ii]; 
             makealanine2(geo,geo.pos,0,nend[ii],phi[ii],psi0[ii]);
             for(k=0;k<=2;k++) // Ca to C vector
               rotv[ii][k] = geo.pos[ nend[ii] ][2][k]-geo.pos[ nend[ii]][1][k];
             mag = sqrt( dot(rotv[ii],rotv[ii]) );  
             for(k=0;k<=2;k++) 
               rotv[ii][k]=rotv[ii][k]/mag;
             for(np[ii]=0;np[ii]<mrot[ii];np[ii]++){
               psi[ii]= psi0[ii]+delrot[ii]*np[ii];
               rotonca(geo,geo.pos,fpos,nend[ii],rotv[ii],geo.pos[nend[ii]][2],psi[ii]); //rotate o,n,and ca about rotv
               tmain = mpt[ii-1][0];
               tpt = mpt[ii-1][1];
               testpt[1] = tpt;
               maincadist(geo,fpos[nend[ii]+1][4],tmain,tpt,nmax,mpt[ii],
                                 &nndist,pos,&bdist);
ccdist = 0;
               maincadist(geo,fpos[nend[ii]+1][1],tmain,tpt,nmax,mpt[ii],
                                 &mdist,pos,&bdist);
#if 1
if(mdist< mtempmin[2]){
   mtempmin[2]= mdist;
}
if(nndist< ntempmin[2]){
   ntempmin[2]= nndist;
}
#endif
               if( (mdist< geo.camaindist)&&(ccdist< geo.camaindist)&&
                   (nndist< geo.camaindist)  )   { //  test another residue
//printf(" 2 successful\n");
//printf(" cadist = %5.3f  ndist = %5.3f  cdist = %5.3f\n",mdist,nndist,ccdist);
                 if(nadd==2){
                   sortbest(fit1,fpos,den1,geo,best,phi,psi,mpt,&iquit,(mdist+nndist)/2.0);
                   if(iquit==1)
                      goto quitrout;
                   goto nextii2;
                 }

                  makealanine2(geo,geo.pos,0,nend[ii],phi[ii],psi[ii]); //if reach here  add good residue
                  ii=2; // third additon 
                  phi0[ii]=psi0[ii] = 0.0; // set phi psi = 0
                  for(nh[ii] =0;nh[ii]<mrot[ii];nh[ii]++){// rotate phi in steps of 10 degrees - find max C den and corresponding Cb den
//printf("nh[%d] = %d\n",ii,nh[ii]);
                     phi[ii]= phi0[ii]+delrot[ii]*nh[ii]; 
                     makealanine2(geo,geo.pos,0,nend[ii],phi[ii],psi0[ii]);
                     for(k=0;k<=2;k++) // Ca to C vector
                       rotv[ii][k] = geo.pos[ nend[ii] ][2][k]-geo.pos[ nend[ii]][1][k];
                     mag = sqrt( dot(rotv[ii],rotv[ii]) );  
                     for(k=0;k<=2;k++) 
                       rotv[ii][k]=rotv[ii][k]/mag;
                     for(np[ii]=0;np[ii]<mrot[ii];np[ii]++){
//printf("np[%d] = %d\n",ii,np[ii]);
                        psi[ii]= psi0[ii]+delrot[ii]*np[ii];
                     rotonca(geo,geo.pos,fpos,nend[ii],rotv[ii],geo.pos[nend[ii]][2],psi[ii]); //rotate o,n,and ca about rotv
                        tmain = mpt[ii-1][0];
                        tpt = mpt[ii-1][1];
                        testpt[2] = tpt;
                        maincadist(geo,fpos[nend[ii]+1][4],tmain,tpt,nmax,mpt[ii],
                                 &nndist,pos,&bdist);
ccdist = 0.0; 
                         maincadist(geo,fpos[nend[ii]+1][1],tmain,tpt,nmax,mpt[ii],
                                 &mdist,pos,&bdist);
if(mdist< mtempmin[3]){
   mtempmin[3]= mdist;
}
if(nndist< ntempmin[3]){
   ntempmin[3]= nndist;
}
                         if( (mdist< geo.camaindist)&&(ccdist< geo.camaindist)&&
                              (nndist< geo.camaindist)  )   { //  test another residue
                         if(nadd==3){
                           sortbest(fit1,fpos,den1,geo,best,phi,psi,mpt,
                                     &iquit,(mdist+nndist)/2.0);
                           if(iquit==1)
                              goto quitrout;
                           goto nextii3;
                         }

                         makealanine2(geo,geo.pos,0,nend[ii],phi[ii],psi[ii]); //if reach here  add good residue
                         ii=3; // third additon 
                         phi0[ii]=psi0[ii] = 0.0; // set phi psi = 0
                         for(nh[ii] =0;nh[ii]<mrot[ii];nh[ii]++){// rotate phi in steps of 10 degrees - find max C den and corresponding Cb den
//printf("nh[%d] = %d\n",ii,nh[ii]);
                            phi[ii]= phi0[ii]+delrot[ii]*nh[ii]; 
                            makealanine2(geo,geo.pos,0,nend[ii],phi[ii],psi0[ii]);
                            for(k=0;k<=2;k++) // Ca to C vector
                               rotv[ii][k] = geo.pos[ nend[ii] ][2][k]-geo.pos[ nend[ii]][1][k];
                            mag = sqrt( dot(rotv[ii],rotv[ii]) );  
                            for(k=0;k<=2;k++) 
                               rotv[ii][k]=rotv[ii][k]/mag;
                            for(np[ii]=0;np[ii]<mrot[ii];np[ii]++){
//printf("np[%d] = %d\n",ii,np[ii]);
                              psi[ii]= psi0[ii]+delrot[ii]*np[ii];
                              rotonca(geo,geo.pos,fpos,nend[ii],rotv[ii],geo.pos[nend[ii]][2],psi[ii]); //rotate o,n,and ca about rotv
                              maincadist(geo,fpos[nend[ii]+1][1],mpt[ii-1][0],mpt[ii-1][1],nmax,mpt[ii],
                                   &mdist,pos,&bdist);
                              if(mdist< geo.camaindist){ //  test another residue or sort and add to set
                               if(nadd==4){
                                 sortbest(fit1,fpos,den1,geo,best,phi,psi,
                                       mpt,&iquit,(mdist+nndist)/2.0);
                                 if(iquit==1)
                                    goto quitrout;
                               }
                             } // mdist
                          }// np[ii]
                        }// nh[ii]
                         ii=2;
                     } // mdist
                     nextii3:; // jump to here if just add 3 residue each step
                     }// np[ii]
                   }// nh[ii]
                   ii=1;
                 } // mdist
                 nextii2:; // jump to here if just add 2 residue each step
              }// np[ii]
            }// nh[ii]
            ii=0;
         } // mdist
         nextii1:; // jump to here if just add 1 residue each step
      } // np[ii]
    } // nh[ii]
    quitrout:;
//    if(verbose>=2){
//    add0types(best);
    if(verbose>=1){
      printf("SORTED TYPES: num = %d: \n",best[0].main[0]);
      for(ii=1;ii<=best[0].main[0];ii++){//each stored type 
         printf("%d: type:",ii);
         for(i=0;i<nadd;i++) // each addition         
             printf("  %d",best[ii].type[i]);
         printf(" qual = %5.3f nbad = %d\n",best[ii].qual,best[ii].nbad);
      }//ii
    }
time1 = clock();
if(verbose>=2)
  printf(" CPU time = %d\n",time1-time0);
printf(" minmdist = %5.3f %5.3f %5.3f  minndist= %5.3f %5.3f %5.3f  \n",mtempmin[1],mtempmin[2],mtempmin[3],
      ntempmin[1],ntempmin[2],ntempmin[3]);
   geo.camaindist = origmaindist;  // set to orig. value when leave
}
#endif


void refineandrank(struct geometry& geo,struct griddata& den1,
      struct pepfit fit1,struct besttype best[NADDTYPE],
      int *nrank,int val[2],struct pepfit fits2[],
      int testpt[],int nmax)
/*
NEW - returns ranking for best density for each residue added
New uses struct best[]
  Only fits with qual >0 are included in best[]
  Uses the good fits found by sortn.
  For each constaint type - if qual>0, 
   refine with tordyn, and if good enough,
   add to ranked sets
*/ 
{
   int i,ii,nbad;
   int useset = 1; // means use set when refining
   int seq0,restype,bestrotn,resn,bestnbad,bestbadset,end1;
   int mpt[3],nbadside,totside;
   float phi,psi,lastden,aveden,avecb2,avedenside;
   float firstden;
   float bestqual;
   float mdist,qualmain,qualside,fbadtot,pos[3],fbadside,bdist;
   int countbad,count;

   seq0 = -fit1.val1;
   if(tordata.useboneconstr ==1){ // If add boneconstrint in dynamics
     for(i=0;i<5;i++)
        tordata.testpt[i] = testpt[i];
     tordata.nmax = nmax;
     tordata.nadd = nadd;
   }
//   countbad = 4; // the number end atoms used in avedenregnadd
//   countbad = 2*nadd; // the number end atoms used in avedenregnadd
   countbad = geo.torend; // the number end atoms used in avedenregnadd
   allocatefit(fitlim,fit1.nres +nadd+2); // make large enough so can hold new fit
   copyfittofit(fit1,fitlim,0,fit1.nres);
   for(ii=1;ii<=best[0].main[0];ii++){ // each constraint type
        if(verbose>=1)
           printf("REFINING CONSTRAINT TYPE %d (total = %d)\n",ii,best[0].main[0]);
        fitend = fit1.nres-1;
        fitbeg = fitend - geo.endmove -4;
        if(fitbeg<0)
           fitbeg = 0;
        geo.main[0] = fit1.main[0];
        geo.main[1] = fit1.main[1];
        geo.main[2] = fit1.main[2];
        copyregfittogeo(geo,fit1,fitbeg,fitend); //NEW
        end1 = geo.nend[1];
        for(i=0;i<nadd;i++){ // each addition
           phi = best[ii].phi[i];
           psi = best[ii].psi[i];
           makealanine2(geo,geo.pos,0,end1+i,phi,psi); 
           restype = geo.seq[ seq0 + fit1.nres-1+i-geo.fnum];

           geo.nend[1]= end1+i+1;
           resn = geo.nend[1]-1; // the side chain that is being built
           if(restype==GLY)
              geo.rtype[geo.nend[1]-1]= 0; 
           else
              geo.rtype[geo.nend[1]-1] = best[ii].type[i];
           bestrotn = -1; // This means try all rotomers.
           refineonerotgeo(useset,den1,geo,rot,resn,restype,&bestqual,
             &bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
           if(verbose>=1)
               printf("For ii=%d  i=%d resn= %d  rottype = %d    bestnbad = %d bestqual = %5.3f\n",
                ii,i,resn,restype,bestnbad,bestqual);
        }//i
        if( (geo.nend[1] - countbad)<geo.nend[0]+1){ //  want to be able to find aveden and nbad over last 2*nadd residues 
           printf("ERROR: geo.nend[1] - countbad)<geo.nend[0]+1\n");
           exit(0);
        }
        geo.nend[1] = end1+nadd;
        tordynbat(30,den1,geo,tordata,&aveden,&nbad);// Run tordyn- BE CAREFUL quitcond must be >5 if use detalilden below 
        maincadist(geo,geo.pos[geo.nend[1]][1],geo.main[0],geo.main[1],
               nadd*nmax,mpt,&mdist,pos,&bdist); // distance from Ca to main
        for(i=0;i<3;i++)
           geo.main[i]=mpt[i];
        avedenregcb(den1,geo,geo.pos,&aveden,&avecb2,&nbad,
                  geo.nend[1]-nadd,geo.nend[1]); // Find average MAINCHAIN + Cb density for last nadd residues
        if(verbose>=1)
              printf("WITH DYN: TYPE = %d %d %d  lastden = %5.3f nbad = %d npoor = %d aveco=%5.3f\n",
                       best[ii].type[0],best[ii].type[1],best[ii].type[2],lastden,nbad,geo.npoor,avecb2);
        qualmain= aveden - nbad*maidfile.denscale;
        qualside =0.0;
        fbadtot = nbad;
        count = 0; // counts residues with sides
        for(i=0;i<nadd;i++){
          if(restype !=GLY){
             count = count+1;
             detaildenoneside(den1,geo,tordata,tordata.nunit1-i,tordata.nunit1-i,&avedenside,&nbadside,
                &firstden,&lastden,&totside); // finddensity for last(=geo.endmove) side chain - same call used in resqual
             if(totside>=4)
                 fbadside = nbadside*4.0/totside;
             else
                 fbadside = nbadside;
             qualside = qualside +avedenside - fbadside*maidfile.denscale;
             fbadtot = fbadtot + fbadside;
         }
       }
       if(count ==0)
           qualside = -10000.0; // indicator that all sides are glycine
       else
           qualside = qualside/count;
       if(verbose>=1)
         printf("After tordynbat qualmain = %5.3f  qualside = %5.3f fbadtot=%5.3f mdist = %5.3f\n",
           qualmain,qualside,fbadtot,mdist);
       if(*nrank< NSAVEFIT-1){// number of saved fits
             *nrank = *nrank+1;
       }
       rankfitside2(nrank,fits2,geo,qualmain,qualside,fbadtot,mdist);
  }// ii
}


void multcbone(struct griddata& den1,struct pepfit fits[NSAVEFIT],
         struct geometry& geo,int *totadd,int iconnect,
         int *lastrank)
/*
On input: fits[1,2,.. lastrank]  are the the saved fits that are being
   extended.  iconnect is the total # of residues needed to reach metend.
On output:  new set of fits[1..]; totadd is the total number of residues
    added from the start of initial extension.
 
 New version - Trys adding nadd residues before refinement:
  1) For each addition - Try all angles and just use those for which 
      end Ca is within geo.camaindist from a main bone,
  2) Keep track of angles for each addition and edit the set to reduce 
       successful sets.
  3) Run tordyn.
*/
{
    int i,j,imain,pt0,nrot,nrank;
    int nmax,val[2],testpt[5];
    struct besttype best[NADDTYPE];
//    float avecb2,lastden,origrfor,qual,aveside;
   struct pepfit fits2[NSAVEFIT];
   char rankfile[10],sfit[10],ranknum[10];

   for(i=1;i<NSAVEFIT;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);
   }

    fitend = fits[1].nres-1; // all fits in fits should be same length
    fitbeg = fitend - geo.endmove -3;
    if(fitbeg <0)
       fitbeg = 0;
    nrot = 36;
//    nrot = 50;
    nmax = geo.gnmax; // number of bone pts to move out to find closept for each residue added
    nrank = 0; // number of saved fits for this extension
    for(i=1;i<= *lastrank;i++){
       if(verbose>=1)
          printf("WORKING ON RANKED FIT %d (total = %d) nadd = %d\n",i,*lastrank,nadd);
       copyregfittogeo(geo,fits[i],fitbeg,fitend); //NEW
       for(j=0;j<3;j++)
           geo.main[j]=fits[i].main[j];
       imain = geo.main[0]; // main that is being fit
       pt0 = geo.main[1]; // last pos on main
//       if(verbose>=2)
          printf("Initial imain = %d  ipt = %d\n",imain,pt0);
       addnresidue(den1,geo,fits[i],imain,pt0,nmax,nrot,best,testpt); // add n residues 
       refineandrank(geo,den1,fits[i],best,&nrank,val,fits2,testpt,nmax);
    }
    if(nrank>0){ // Then successfully added one  set of nadd residue
         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
         *totadd = *totadd + nadd;
   } //nrank>0 -  savefit
   else
      *lastrank = 0;
  for(i=1;i<NSAVEFIT;i++) 
      freepepfit(fits2[i]);
}
