/* assignseq.c++  Routines for assigning sequence to prevously fit segments.
9dec99:  Problem - when use geo.set  - must check to see if atom out of range.
      Also - should skip over gly and pro  - not loop through stuff
27oct99:  Major Revision:
    Use geo.set (set useset = 1)  so that can test for bad contacts.
    Big improvement for  Arg and Lys 

30aug99:  Am debugging routine avedengeooneside:  
              Works OK but counts some atmos twice.
6aug99: Modified so that program does not copy entire fit to geo
          but instead just copies one fit res to geo.
       Seems to give identical result - (BUT I AM NOT SURE)
         If questions - see backup assignseq.c5aug99
5aug99:  Program seems to work well
4aug99
  Works OK - after much debugging
  However, in fit pcab:  Q 130   -  fum 163  - there is an
      arginine that overlaps another res.
      TRY TO FIX!
29jul99 - working on being able to save and output the fits that have ahd
   residue sidechains assigned.
  Procedure:
    Keep modifying current fit (and saved fit)
    When savefit, automatically write backup fit, this allows one to 
       readbackupfit and start over if do not like.
Some rules for assigning sequences (2dec98):
  Using Powell Minimization seems to work better and faster than torsion dynamics.
    (although nothing seems to be wrong with the dynamics routine)
  1) Using geo.minrigid = 0.5*denscale and 
       qual = (bestden) - 0.5*(nbad)*maidfile.denscale - 0.5*(badset)*maidfile.denscale;
     Seems to provide good discrimination.
  2) Value of aveden >40 seems to be a good fit.
  3) For residues like PHE (TYR,HIS) - if aveden>40, the stonng indicator.
  4) PHE can have good fit with TRP, but not vice -versa,
       Thus compare TRP and PHE, if TRP about = to PHE, implies
      it is a PHE
  5) For some reason, long chains (ARG) can get pretty good fit in
       lots of places  - cannot really descriminate.
     Do not try to determine if ARG
  6) MET : Use criteria that den atom 7 > 2.0(?)*atom 6
     Example:
        For ifit = 3 resn= 5  rottype = 12   Best rot = 4: bestnbad = 0 bestqual = 38.938
        atomden 6,7,8 = 23.070 68.695 36.218

*/
#include "pdbtolev.h"
#include "drawwindow.h"  // needed only for maidfile.denscale
#include "math.h"
#include "trace.h" //NOTE  fit and fitnum are defined globally in trace.h; includes a number of "extern .... " lines
#include "peptide.h"
#include "tordyn.h"
#include "extendfit.h"
#include "initialize.h"
#include "editbones.h"
#include <string.h>
#include "assignseq.h"
#include "fittopdb.h"
#include "rankside.h"
#include "symmate.h"


extern int verbose; // if = 1, print a lot of stuff
extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version
extern int side[AA][NBRN][NPTS];  // defined globally (originally in pdbtolev.c++)  because needed in assignseg.c++ 
extern int brdraw[AA][NBRN]; // defined globally because MAYBE?? needed in assignseg.c++

extern struct maiddatafile maidfile; // needed only for maidfile.denscale
extern struct tordyn tordata;
extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern int fitnum[NFIT];// see trace.c++ for definition
extern struct extbones ext; // All the data for edited bones used to extend one end of a given fit
extern struct rotomer rot[30];// rotomer data for each amino acid
extern char*  nameres[30];
extern int data[AA+4][50]; 

void findside_maindisttor(struct geometry& geo,
       int ires,int iside,float *mindist);
   
void removefit(struct pepfit fit[],int ifit)
/*
   Removes ifit from set of fits and writes to fit file
*/
{
   int i;

   for(i= ifit+1;i<=fitnum[0];i++)
      copyfittofit(fit[i],fit[i-1],0,fit[i].nres); // copy higher number fit to previous number
   fitnum[0] = fitnum[0]-1;
   savefit();
}



void removenunordfits(struct pepfit fit[],int del[])
/* Used if fits are not ordered
  First calls orderdel, then removenfits
*/
{
   int order[MAXDEL];
  
   if(del[0]>MAXDEL-2){
     printf("ERROR: del[0]=%d>MAXDEL: setting del[0]=MAXDEL-2\n",del[0]);
     printf("     INCREASE MAXDEL AND RECOMPILE\n");
     del[0]=MAXDEL;
   }
   orderlist(del,order);
   removenfits(fit,order);
}
   
void removenfits(struct pepfit fit[],int del[])
/*
   Removes del[0] fits (del[1,2,..] from set of fits and writes to fit file
  NOTE: IT IS ASSUMED THAT del[1,2,..] is ordered in order of increasing fit #
*/
{
   int i,irem,count;

  
printf("REMOVING %d Fits:",del[0]);
   for(i=1;i<=del[0];i++)
       printf("  %d",del[i]);
   printf("\n");
   del[ del[0]+1] = fitnum[0]+1; // artificial value so routine works
   count = 0; // counts number removed
   for(irem = 1;irem<=del[0];irem++){
      count = count+1;
      for(i= del[irem]+1;i<del[irem+1];i++)
          copyfittofit(fit[i],fit[i-count],0,fit[i].nres); // copy higher number fit to previous number
   }
   fitnum[0] = fitnum[0]-count;
   savefit();
}

void copyonerottofit(float rotpos[NATOMS][3],struct pepfit& fit,int fitloc,int rottype)// copy rotomer side chain to fit
// use count to convert from atomname to side atoms numbers 6,7,8  etc.
{
   int ii,j,k,atomname,count;

   count = 0;
   fit.nside[fitloc] = -rottype;  // minus sign (< 0) indicates have residue side chain, not just series od C atoms
   for(ii=1;ii<=side[rottype][0][0]+1;ii++){ // copy each side chain branch
      for(j=0;j<side[rottype][ii][0];j++){
         atomname=side[rottype][ii][j+1];
         if(atomname>5){ // not part of mainchain or Cb
           count = count+1;
//printf("atomname = %d\n",atomname);
           for(k=0;k<=2;k++)
              fit.pos[fitloc][count+5][k]=rotpos[atomname][k];
         }
         else{
           for(k=0;k<=2;k++)
              fit.pos[fitloc][atomname][k]=rotpos[atomname][k];
         }
      }
   } //ii
}

void copyonerottogeo(float rotpos[NATOMS][3],struct geometry& geo,int geoloc,int rottype)// copy rotomer side chain to fit
// use count to convert from atomname to side atoms numbers 6,7,8  etc.
{
   int ii,j,k,atomname,count;

   count = 0;
   geo.iside[geoloc] = -rottype;  // minus sign (< 0) indicates have residue side chain, not just series od C atoms
   for(ii=1;ii<=side[rottype][0][0]+1;ii++){ // copy each side chain branch
      for(j=0;j<side[rottype][ii][0];j++){
         atomname=side[rottype][ii][j+1];
         if(atomname>5){ // not part of mainchain or Cb
           count = count+1;
//printf("atomname = %d\n",atomname);
           for(k=0;k<=2;k++)
              geo.pos[geoloc][count+5][k]=rotpos[atomname][k];
         }
         else{
           for(k=0;k<=2;k++)
              geo.pos[geoloc][atomname][k]=rotpos[atomname][k];
         }
      }
   } //ii
   if( rottype==GLY)
      geo.resside[geoloc] =4;
   else if( rottype==ALA)
      geo.resside[geoloc] =5;
   else 
      geo.resside[geoloc] =5+count;
}

void rotlib(int aanum,int *rotnum,float chiang[11][3])
/*
New - 13sep99  Include 3 rotomers for PRO
 return # of rotomers (rotnum) and the chiangles  (in Radians)
    for res type aanum)
   listed in order of frequency found
 */
{
   int i,j;

   float conv=PI/180.0;
   switch(aanum){
     case(ALA): /* the first three cases are added just to get correct rotnum*/
        *rotnum=1;
        break;
     case(GLY):
        *rotnum=1;
        break;
#if 0 // These are values in library - Do not use because adds unnecessary complication
     case(PRO):
        *rotnum=3;
        chiang[1][0]= 26.9; /* chi1*/
        chiang[1][1]= -29.4; /* chi2*/
        chiang[2][0]= -21.8;
        chiang[2][1]= 31.2;
        chiang[3][0]= 0.3;
        chiang[3][1]= -0.8;
        break;
#endif
#if 1 // Just use chiang = 0.0 for chiang [1] and [2]
     case(PRO):
        *rotnum=1;
        break;
#endif
     case(CYS):
        *rotnum=3;
        chiang[1][0]= -65.; /* chi1*/
        chiang[2][0]= -178.;
        chiang[3][0]= 64.;
        break;
     case(SER):
        *rotnum=3;
        chiang[1][0]= 65.; /* chi1*/
        chiang[2][0]= -63.;
        chiang[3][0]= -178.;
        break;
     case(THR):
        *rotnum=3;
        chiang[1][0]= -61.; /* chi1*/
        chiang[2][0]= 65.;
        chiang[3][0]= -160.;
        break;
     case(VAL):
        *rotnum=3;
        chiang[1][0]= 173.; /* chi1*/
        chiang[2][0]= -62.;
        chiang[3][0]= 72.;
        break;
     case(ASN): // interchange 4 and 3
        *rotnum=4;
        *rotnum=3; //Number that are sig. diff
        chiang[1][0]= -72.; /* chi1*/
        chiang[1][1]= -44; /* chi2*/
        chiang[2][0]= -171.;
        chiang[2][1]= 12;
        chiang[3][0]= 66;
        chiang[3][1]= 12;
        chiang[4][0]= -69.;
        chiang[4][1]= 131;
        break;
     case(ASP):
        *rotnum=3;
        chiang[1][0]= -70.; /* chi1*/
        chiang[1][1]= -32; /* chi2*/
        chiang[2][0]= -171.;
        chiang[2][1]= -171;
        chiang[3][0]= 62.;
        chiang[3][1]= 4;
        break;
     case(HIS):// interchange 5 and 3
        *rotnum=6;
        *rotnum=3; // Only include rotomers that are significantly different
        chiang[1][0]= -62.; /* chi1*/
        chiang[1][1]= -84; /* chi2*/
        chiang[2][0]= -177.;
        chiang[2][1]= 68;
        chiang[3][0]= 64; 
        chiang[3][1]= -95;
        chiang[4][0]= -168;
        chiang[4][1]= -108;
        chiang[5][0]= -67.;
        chiang[5][1]= 97;
        chiang[6][0]=68;
        chiang[6][1]=81;
        break;
#if 0 // PROBLEM = ILE DOES NOT SATISFY ASSUMPTION IN TORDYN THAT ONLY LAST RESIDUE HAS MORE THAN ONE ATOM 
     case(ILE): // OLd version - used to 8jun00
        *rotnum=5;
        chiang[1][0]= -62.; /* chi1*/
        chiang[1][1]= -163; /* chi2*/  (= +197)
        chiang[2][0]= 63.;
        chiang[2][1]= 164;
        chiang[3][0]= -57;
        chiang[3][1]= -59;
        chiang[4][0]= -168;
        chiang[4][1]= 161;
        chiang[5][0]= -156; 
        chiang[5][1]= 72;
        break;
/* Probability (3rd column, chi1 and chi2 in last 4 columns) of different rotomers (from Dunbrack)
 4 highest probabilities are:
ILE 1 2 0 0   1399   1279   12.59   0.27   91.34   0.61   61.5  8.6    170.6 11.3  
ILE 2 2 0 0   1008    678    6.68   0.20   67.19   1.21 -171.5 12.2    167.7 12.3  
ILE 3 2 0 0   7751   5861   57.67   0.40   75.61   0.40  -63.9  7.7    169.9 10.7  
ILE 3 3 0 0   7751   1628   16.02   0.30   21.01   0.38  -59.2  8.6    -61.2 11.7
 75%  have second chi = 2   = 168 
       r1,r2,r3,r4     Conformation       chi range
        -----------     ------------    -----------------
        1               g+                 0 <= chi < 120
        2               t                120 <= chi < 240
        3               g-              -120 <= chi <   0     
ILE 1 1 0 0   1399    105    1.04   0.08    7.54   0.58   62.0 26.2     91.9 15.0  
ILE 1 2 0 0   1399   1279   12.59   0.27   91.34   0.61   61.5  8.6    170.6 11.3  
ILE 1 3 0 0   1399     15    0.15   0.03    1.12   0.23   82.4 17.7    -82.1 18.8  
ILE 2 1 0 0   1008    313    3.09   0.14   31.06   1.19 -167.2 15.7     68.7 13.9  
ILE 2 2 0 0   1008    678    6.68   0.20   67.19   1.21 -171.5 12.2    167.7 12.3  
ILE 2 3 0 0   1008     17    0.17   0.03    1.75   0.34 -163.3 26.6    -80.3 25.8  
ILE 3 1 0 0   7751    262    2.58   0.13    3.39   0.17  -72.8 18.9     81.0 22.6  
ILE 3 2 0 0   7751   5861   57.67   0.40   75.61   0.40  -63.9  7.7    169.9 10.7  
ILE 3 3 0 0   7751   1628   16.02   0.30   21.01   0.38  -59.2  8.6    -61.2 11.7
*/
#endif
#if 1 // // NEW version - used after 8jun00 - assume chiang[2] is always = 164
    case(ILE):
        *rotnum=3;
        chiang[1][0]= -63.9; /* chi1*/ 
        chiang[2][0]= 61.5;
        chiang[3][0]= -171.5;  
        for(i=1;i<=3;i++) // This is not actually used in buildrot()
           chiang[i][1] = 169.0;
        break;
#endif
       case(LEU):
        *rotnum=6;
        chiang[1][0]= -62.; /* chi1*/
        chiang[1][1]= 170; /* chi2*/
        chiang[2][0]= -176;
        chiang[2][1]= 66;
        chiang[3][0]= -99;
        chiang[3][1]= 25;
        chiang[4][0]= -168;
        chiang[4][1]= -161;
        chiang[5][0]= -83; 
        chiang[5][1]= -80;
        chiang[6][0]=61;
        chiang[6][1]=85;
        break;
     case(AA3): // artificial 3 atom straight - same as LEU
        *rotnum=6;
        chiang[1][0]= -62.; /* chi1*/
        chiang[1][1]= 170; /* chi2*/
        chiang[2][0]= -176;
        chiang[2][1]= 66;
        chiang[3][0]= -99;
        chiang[3][1]= 25;
        chiang[4][0]= -168;
        chiang[4][1]= -161;
        chiang[5][0]= -83; 
        chiang[5][1]= -80;
        chiang[6][0]=61;
        chiang[6][1]=85;
        break;
     case(PHE):
        *rotnum=3;
        chiang[1][0]= 66; /* chi1*/
        chiang[1][1]= 89; /* chi2*/
        chiang[2][0]= -177;
        chiang[2][1]= 73;
        chiang[3][0]= -66;
        chiang[3][1]= 102;
        break;
     case(TRP):
        *rotnum=7;
        chiang[1][0]= -178; chiang[1][1]= 82;         
        chiang[2][0]= -67;chiang[2][1]= 85;
        chiang[3][0]= -174;chiang[3][1]= -101;        
        chiang[4][0]= 64;chiang[4][1]= 93;        
        chiang[5][0]= 57; chiang[5][1]= -86;        
        chiang[6][0]= -74;chiang[6][1]= -86;
        chiang[7][0]= -73;chiang[7][1]= 4;
     case(TYR):
        *rotnum=6;
        chiang[1][0]= -64; chiang[1][1]= 102;         
        chiang[2][0]= -66;chiang[2][1]= -87;
        chiang[3][0]= -178;chiang[3][1]= 74;        
        chiang[4][0]= -177;chiang[4][1]= -103;        
        chiang[5][0]= 65; chiang[5][1]= -96;        
        chiang[6][0]= 62;chiang[6][1]=89;        
        break;
     case(GLN)://Interchange 2 and 9
        *rotnum=10;
        *rotnum=5; //number significantly different
        chiang[1][0]= -63; chiang[1][1]= 178;chiang[1][2]= -3;         
        chiang[2][0]= 69;chiang[2][1]= -173;chiang[2][2]= -10; 
        chiang[3][0]= -68;chiang[3][1]= -64; chiang[3][2]= -42;        
        chiang[4][0]= -173;chiang[4][1]= 180;chiang[4][2]= 1;         
        chiang[5][0]= -172; chiang[5][1]= 68;chiang[5][2]= 32;         
        chiang[6][0]= -66;chiang[6][1]= -70;chiang[6][2]= 140; 
        chiang[7][0]= -174;chiang[7][1]= -177 ;chiang[7][2]= 179; 
        chiang[8][0]= -170; chiang[8][1]= 68;chiang[8][2]= -124;         
        chiang[9][0]= -68;chiang[9][1]= -179;chiang[9][2]= -176; 
        chiang[10][0]= 64;chiang[10][1]= -169;chiang[10][2]= 116; 
        break;
     case(GLU):// Interchange 8 and 2, 7 and 5 , 9 and 6
        *rotnum=10;
        *rotnum=6;//number significantly different
        chiang[1][0]= -66; chiang[1][1]= 178;chiang[1][2]= 174;         
        chiang[2][0]= 72; chiang[2][1]= 178;chiang[2][2]= 41;         
        chiang[3][0]= -64;chiang[3][1]= -65; chiang[3][2]= -45;        
        chiang[4][0]= -174;chiang[4][1]= 176;chiang[4][2]= -173;         
        chiang[5][0]= -65;chiang[5][1]= -79 ;chiang[5][2]= 8; 
        chiang[6][0]= 50;chiang[6][1]= -87;chiang[6][2]= -161; 
        chiang[7][0]= 179; chiang[7][1]= 175;chiang[7][2]= 8;         
        chiang[8][0]= -69;chiang[8][1]= 180;chiang[8][2]= -15; 
        chiang[9][0]= -175;chiang[9][1]= 70;chiang[9][2]= 57; 
        chiang[10][0]= 65;chiang[10][1]= 176;chiang[10][2]= -120; 
        break;
     case(AA4)://artficial straint 4 atom - same as GLU 
        *rotnum=10;
        *rotnum=6;//number significantly different
        chiang[1][0]= -66; chiang[1][1]= 178;chiang[1][2]= 174;         
        chiang[2][0]= 72; chiang[2][1]= 178;chiang[2][2]= 41;         
        chiang[3][0]= -64;chiang[3][1]= -65; chiang[3][2]= -45;        
        chiang[4][0]= -174;chiang[4][1]= 176;chiang[4][2]= -173;         
        chiang[5][0]= -65;chiang[5][1]= -79 ;chiang[5][2]= 8; 
        chiang[6][0]= 50;chiang[6][1]= -87;chiang[6][2]= -161; 
        chiang[7][0]= 179; chiang[7][1]= 175;chiang[7][2]= 8;         
        chiang[8][0]= -69;chiang[8][1]= 180;chiang[8][2]= -15; 
        chiang[9][0]= -175;chiang[9][1]= 70;chiang[9][2]= 57; 
        chiang[10][0]= 65;chiang[10][1]= 176;chiang[10][2]= -120; 
        break;
     case(MET):
        *rotnum=10;
        chiang[1][0]= -66; chiang[1][1]= -67;chiang[1][2]= -67;         
        chiang[2][0]= -68;chiang[2][1]= 180;chiang[2][2]= 74; 
        chiang[3][0]= -72;chiang[3][1]= 179; chiang[3][2]= 180;        
        chiang[4][0]= 179;chiang[4][1]= 178;chiang[4][2]= 77;         
        chiang[5][0]= -169; chiang[5][1]= -169;chiang[5][2]= 101;         
        chiang[6][0]= -61;chiang[6][1]= -67;chiang[6][2]= 131; 
        chiang[7][0]= -170;chiang[7][1]= 65;chiang[7][2]= 76; 
        chiang[8][0]= 80; chiang[8][1]= -172;chiang[8][2]= 180;         
        chiang[9][0]= -74;chiang[9][1]= -173;chiang[9][2]= -82; 
        chiang[10][0]= -66;chiang[10][1]= -67;chiang[10][2]= 140; 
        break;
     case(LYS):
        *rotnum=10;
        chiang[1][0]= -70; chiang[1][1]= -179;chiang[1][2]= 180;         
        chiang[2][0]= -172;chiang[2][1]= 180;chiang[2][2]= 180; 
        chiang[3][0]= -61;chiang[3][1]= -71; chiang[3][2]= -178;        
        chiang[4][0]= -173;chiang[4][1]= -174;chiang[4][2]= 77;         
        chiang[5][0]= -174; chiang[5][1]= 70;chiang[5][2]= -174;         
        chiang[6][0]= 70;chiang[6][1]= -172;chiang[6][2]= 179; 
        chiang[7][0]= -69;chiang[7][1]= 178;chiang[7][2]= 79; 
        chiang[8][0]= -71; chiang[8][1]= -172;chiang[8][2]= -72;         
        chiang[9][0]= -170;chiang[9][1]= -176;chiang[9][2]= -74; 
        chiang[10][0]= -71;chiang[10][1]= 90;chiang[10][2]= -168; 
        break;
     case(ARG):
        *rotnum=10;
        chiang[1][0]= -176; chiang[1][1]= 156;chiang[1][2]= -175;         
        chiang[2][0]= -67;chiang[2][1]= -176;chiang[2][2]= -76; 
        chiang[3][0]= -178;chiang[3][1]= 171; chiang[3][2]= 62;        
        chiang[4][0]= -72;chiang[4][1]= 166;chiang[4][2]= 72;         
        chiang[5][0]= -60; chiang[5][1]= -60;chiang[5][2]= 171;         
        chiang[6][0]= 78;chiang[6][1]= 171;chiang[6][2]= 172; 
        chiang[7][0]= -172;chiang[7][1]= -170;chiang[7][2]= -66; 
        chiang[8][0]= -69; chiang[8][1]= 175;chiang[8][2]= 180;         
        chiang[9][0]= -59;chiang[9][1]= -74;chiang[9][2]= -75; 
        chiang[10][0]= 66;chiang[10][1]= 158;chiang[10][2]= 74; 
        break;
    }
    for(i=1;i<= *rotnum;i++) /* convert to radians*/
       for(j=0;j<3;j++)
          chiang[i][j]=conv*chiang[i][j];
}    

void cnext(struct geometry& geo,float c1[3],float c2[3],float c3[3],float al1,
      float al2,float lbond,float chi,float c4[3]) 
/* finds the next atom (C4) in sidechain given preceding 3 atoms (c1-c3).
    the two previous angles (al1,al2), the next bond length (lbond),
   and the chi angle 
Example: C1=N, c2= Ca, C3= Cb, and C4= Cgamma
  al1 = c1-c2-c3 angle, al2= c2-c3-c4 angle
NOTE: al2 has sign: + if c2-c1 and c2-c3 point in same direction, - if opposite
   In this case always use + sign because it is what want for chi=0; then rotate
 (The negative sign would correspond to chi= -180)
Procedure:
   1) First place c4 in same place as c1-c3.
   2) then rotate c4 by chi about c2-c3 ( right hand rule )
Note:  al1 and al2 are the angles when vectors are placed tail to end.
   the angle in routine planevec are when vectors have same origin.
*/
{
   int i;
   float u1[3],u2[3],mag1,mag2,theta,ang;
   float costh,sinth;

   for(i=0;i<=2;i++){
     u1[i]=c2[i]-c1[i];
     u2[i]=c3[i]-c2[i];
   }
   mag1=sqrt(dot(u1,u1));
   mag2=sqrt(dot(u2,u2));
   for(i=0;i<=2;i++){
     u1[i]=u1[i]/mag1;
     u2[i]=u2[i]/mag2;
   }
   theta=PI-al1; 
   costh=cos(theta);
   sinth=1-sqr(costh); /* sin squared */ 
   planevec(u1,u2,costh,sinth,lbond,theta,PI-al2+theta,c3,c4); /*  this
      returns vector =c4-c3 (lenght lbond) in same plane as u1,u2*/
   rotatept(geo,u2,c3,chi,c4); 
/* check */
for(i=0;i<=2;i++)
  u1[i]=c4[i]-c3[i];
ang =acos(dot(u2,u1)/lbond);
mag1=sqrt(dot(u1,u1));
ang=PI-ang;
}

void mnext(float c1[3],float c2[3],float c3[3],float al1,
      float al2,float lbond,float c4[3]) 
/* Given c1-c2-c3, finds atom (c4) that has bond vector c2-c4 that is in 
    sameplane as c1-c2-c3.
al1 = c1-c2-c3 angle
al2=  c1-c2-c4 angle
Note: do not need to rotate by chi angle, since in all cases there is no rotation
Note: use negataive: -(PI-al2)=al2-PI because need negative angle in plane vec for th1
*/
{
   int i;
   float u1[3],u2[3],mag1,mag2,theta;
   float costh,sinth;


   for(i=0;i<=2;i++){
     u1[i]=c2[i]-c1[i];
     u2[i]=c3[i]-c2[i];
   }
   mag1=sqrt(dot(u1,u1));
   mag2=sqrt(dot(u2,u2));
   for(i=0;i<=2;i++){
     u1[i]=u1[i]/mag1;
     u2[i]=u2[i]/mag2;
   }
   theta=PI-al1; 
   costh=cos(theta);
   sinth=1-sqr(costh); /* sin squared */ 
   planevec(u1,u2,costh,sinth,lbond,theta,al2-PI,c2,c4); /*  this
      returns vector =c4-c2 (lenght lbond) in same plane as u1,u2*/
}

  
void buildrot(struct geometry& geo,float n[3],float ca[3],float cb[3],
         int restype,float chiang[3],float rotpos[NATOMS][3])
/*
n and ca and cb are positions of preceding n and Ca atoms
Before Using - must call
Builds the rotopmer
 builds the rotamer for pdb # pdbnum, side chain residue aanum
  For glycine, a new CB is build using ideal angles, etc.
   Note this ideal position use ideal tetraheryl bond (= 109.47) while
           M. Levitt has this bond at 112????
  Note mainchain in old maid is defined as:
    mainchain[ segment ][ pdb ][ 0=N, 1=CA, 2=C ][x,y,z]
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i] and fit.pos
*/
{
    float caton[3],catocb[3],mag1,mag2,al1;
    int k;
    float tang;


    tang=120.0*PI/180.0; /* rot. angle for tetrangonal Cargon*/
    for(k=0;k<=2;k++){
          rotpos[N][k]= n[k];
          rotpos[CA][k]= ca[k];
    }
#if 0 // this was needed in oldmaid because no Cb atom for glycine
    if(strcmp(rname[seg][pdbnum],"GLY")==0){
       idealcb(rotpos[N],rotpos[CA],mainch[seg][pdbnum][2],rotpos[CB]);
    }
    else{
#endif
    for(k=0;k<=2;k++)
          rotpos[CB][k]= cb[k];
    for(k=0;k<=2;k++){ /* same for all rotomers i */
      caton[k]=rotpos[N][k]-rotpos[CA][k]; /* vetor from ca to n */
      catocb[k]=rotpos[CB][k]-rotpos[CA][k]; /* vetor from ca to cb */
    }
//printf("N = %5.3f %5.3f %5.3f\n",rotpos[N][0],rotpos[N][1],rotpos[N][2]);
//printf("Ca = %5.3f %5.3f %5.3f\n",rotpos[CA][0],rotpos[CA][1],rotpos[CA][2]);
//printf("Cb = %5.3f %5.3f %5.3f\n",rotpos[CB][0],rotpos[CB][1],rotpos[CB][2]);
    mag1=sqrt(dot(caton,caton));
    mag2=sqrt(dot(catocb,catocb));
    al1=acos(dot(caton,catocb)/(mag1*mag2));/* the N-Ca-Cb angle */

   switch(restype){
      case(ALA): /* just use the above ideal CB */
         break;
      case(ARG): /* arginine , assume trans (=PI) rotations for last two chi*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.ccc,geo.ccn,geo.cn,chiang[2],rotpos[NE]);
            cnext(geo,rotpos[CG],rotpos[CD],rotpos[NE],geo.ccn,geo.anc,geo.an,PI,rotpos[CZ]);
            cnext(geo,rotpos[CD],rotpos[NE],rotpos[CZ],geo.anc,geo.nan,geo.an,PI,rotpos[NH1]);
            mnext(rotpos[NE],rotpos[CZ],rotpos[NH1],geo.nan,geo.nan,geo.an,rotpos[NH2]);
         break;
      case(ASN):/* asparagine */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.acc,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.cao,geo.ao,chiang[1],rotpos[OD1]);
            mnext(rotpos[CB],rotpos[CG],rotpos[OD1],geo.cao,geo.can,geo.an,rotpos[ND2]);
         break;
      case(ASP):/* aspartate */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.acc,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.cao,geo.ao,chiang[1],rotpos[OD1]);
            mnext(rotpos[CB],rotpos[CG],rotpos[OD1],geo.cao,geo.cao,geo.ao,rotpos[OD2]);
         break;
      case(CYS):/* cysteine */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccs,geo.cs,chiang[0],rotpos[SG]);
         break;
      case(GLU):/* glutamate */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.acc,geo.ac,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.acc,geo.cao,geo.ao,chiang[2],rotpos[OE1]);
            mnext(rotpos[CG],rotpos[CD],rotpos[OE1],geo.cao,geo.cao,geo.ao,rotpos[OE2]);
         break;
      case(AA4):// artificial residue used to sample density
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.acc,geo.ac,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.acc,geo.cao,geo.ao,chiang[2],rotpos[OE1]);
//            mnext(rotpos[CG],rotpos[CD],rotpos[OE1],geo.cao,geo.cao,geo.ao,rotpos[OE2]);
         break;
      case(GLN):/* glutamine */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.acc,geo.ac,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.acc,geo.cao,geo.ao,chiang[2],rotpos[OE1]);
            mnext(rotpos[CG],rotpos[CD],rotpos[OE1],geo.cao,geo.can,geo.an,rotpos[NE2]);
         break;
      case(GLY):
         break;
      case(HIS):/* histadine  Note third chi = PI, rest of ring = 0*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.phe,geo.bc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bcc,geo.cbn,geo.bn,chiang[1],rotpos[ND1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[ND1],geo.cbn,geo.bnb,geo.bn,PI,rotpos[CE1]);
            cnext(geo,rotpos[CG],rotpos[ND1],rotpos[CE1],geo.bnb,geo.nbn,geo.bn,0.0,rotpos[NE2]);
            cnext(geo,rotpos[ND1],rotpos[CE1],rotpos[NE2],geo.nbn,geo.bnb,geo.bn,0.0,rotpos[CD2]);
         break;
      case(ILE):// Major problem - tordyn assumes only terminal group has more than one atom
               //  Thus, treat all 4 terminal atoms as 1 fixed group - ASSUME chiang[i][1] = 169.0;
            chiang[1]= 169.0*PI/180.0;
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG1]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG1],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD1]);
#if 0
            mnext(rotpos[CA],rotpos[CB],rotpos[CG1],geo.ccc,ccc,geo.cc,rotpos[CG2]);
#endif
#if 1
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0]-tang,rotpos[CG2]);
#endif
         break;
      case(LEU):
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1]+tang,rotpos[CD2]);
         break;
      case(AA3): // artificial residue, used to sample density
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD1]);
//            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1]+tang,rotpos[CD2]);
         break;
      case(LYS): /* assume trans rot for last chi*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.ccc,geo.ccc,geo.cc,chiang[2],rotpos[CE]);
            cnext(geo,rotpos[CG],rotpos[CD],rotpos[CE],geo.ccc,geo.ccn,geo.cn,PI,rotpos[NZ]);
         break;
      case(MET): /* assume trans rot for last chi*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccs,geo.cs,chiang[1],rotpos[SD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[SD],geo.ccs,geo.csc,geo.cs,chiang[2],rotpos[CE]);
         break;
      case(PHE): /* use aaa where want aac (not in table( */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.phe,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.aaa,geo.aa,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD1],geo.aaa,geo.aaa,geo.aa,PI,rotpos[CE1]);
            cnext(geo,rotpos[CG],rotpos[CD1],rotpos[CE1],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CZ]);
            cnext(geo,rotpos[CD1],rotpos[CE1],rotpos[CZ],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CE2]);
            cnext(geo,rotpos[CE1],rotpos[CZ],rotpos[CE2],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CD2]);
         break; 
#if 0     
      case(PRO):
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.bbb,geo.bb,chiang[0],rotpos[CG]);
         cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bbb,geo.bbb,geo.bb,chiang[1],rotpos[CD]);
         break;
#endif
#if 1   // New = fix chiang  0   
      case(PRO):
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.bbb,geo.bb,0.0,rotpos[CG]);
         cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bbb,geo.bbb,geo.bb,0.0,rotpos[CD]);
         break;
#endif
      case(SER):
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccv,geo.cv,chiang[0],rotpos[OG]);
         break;
      case(THR):
           cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccv,geo.cv,chiang[0],rotpos[OG1]);
#if 0
           mnext(rotpos[CA],rotpos[CB],rotpos[OG1],geo.ccv,ccv,geo.cv,rotpos[CG2]);
#endif
           cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccv,geo.cv,chiang[0]+tang,rotpos[CG2]);
         break;
      case(TRP): /* Note third chi = PI */
//            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.bcc,geo.bc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.phe,geo.bc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bcc,geo.bbc,geo.bb,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD1],geo.bbc,geo.bbn,geo.bn,PI,rotpos[NE1]);
            cnext(geo,rotpos[CG],rotpos[CD1],rotpos[NE1],geo.bbn,geo.bng,geo.gn,0.0,rotpos[CE2]);
            cnext(geo,rotpos[CD1],rotpos[NE1],rotpos[CE2],geo.bng,geo.ggn,geo.gg,0.0,rotpos[CD2]);
            mnext(rotpos[CE2],rotpos[CD2],rotpos[CG],geo.bgg,geo.agg,geo.ag,rotpos[CE3]);
            cnext(geo,rotpos[CE2],rotpos[CD2],rotpos[CE3],geo.agg,geo.aag,geo.aa,0.0,rotpos[CZ3]);
            cnext(geo,rotpos[CD2],rotpos[CE3],rotpos[CZ3],geo.aag,geo.aaa,geo.aa,0.0,rotpos[CH2]);
            cnext(geo,rotpos[CE3],rotpos[CZ3],rotpos[CH2],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CZ2]);
         break;
      case(TYR): /* use aaa where want aac (not in table) */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.tyr,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.aaa,geo.aa,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD1],geo.aaa,geo.aaa,geo.aa,PI,rotpos[CE1]);
            cnext(geo,rotpos[CG],rotpos[CD1],rotpos[CE1],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CZ]);
            cnext(geo,rotpos[CD1],rotpos[CE1],rotpos[CZ],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CE2]);
            cnext(geo,rotpos[CE1],rotpos[CZ],rotpos[CE2],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CD2]);
            mnext(rotpos[CE1],rotpos[CZ],rotpos[CE2],geo.aaa,geo.aav,geo.av,rotpos[OH]);

         break;      
      case(VAL):
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG1]);
#if 0
         mnext(rotpos[CA],rotpos[CB],rotpos[CG1],geo.ccc,ccc,geo.cv,rotpos[CG2]);
#endif
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0]+tang,rotpos[CG2]);
         break;
            
    }
}

  
void buildrotring(struct geometry& geo,float n[3],float ca[3],float cb[3],
         int restype,float chiang[3],float ringang,float rotpos[NATOMS][3])
/*
NEW - 29jan01 - allows variable ring angle -
  Only called after tordyn 
In assignseq() routines - only use standard buildrot() - use default ringang.
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i] and fit.pos
*/
{
    float caton[3],catocb[3],mag1,mag2,al1;
    int k;
    float tang;


    tang=120.0*PI/180.0; /* rot. angle for tetrangonal Cargon*/
    for(k=0;k<=2;k++){
          rotpos[N][k]= n[k];
          rotpos[CA][k]= ca[k];
    }
#if 0 // this was needed in oldmaid because no Cb atom for glycine
    if(strcmp(rname[seg][pdbnum],"GLY")==0){
       idealcb(rotpos[N],rotpos[CA],mainch[seg][pdbnum][2],rotpos[CB]);
    }
    else{
#endif
    for(k=0;k<=2;k++)
          rotpos[CB][k]= cb[k];
    for(k=0;k<=2;k++){ /* same for all rotomers i */
      caton[k]=rotpos[N][k]-rotpos[CA][k]; /* vetor from ca to n */
      catocb[k]=rotpos[CB][k]-rotpos[CA][k]; /* vetor from ca to cb */
    }
//printf("N = %5.3f %5.3f %5.3f\n",rotpos[N][0],rotpos[N][1],rotpos[N][2]);
//printf("Ca = %5.3f %5.3f %5.3f\n",rotpos[CA][0],rotpos[CA][1],rotpos[CA][2]);
//printf("Cb = %5.3f %5.3f %5.3f\n",rotpos[CB][0],rotpos[CB][1],rotpos[CB][2]);
    mag1=sqrt(dot(caton,caton));
    mag2=sqrt(dot(catocb,catocb));
    al1=acos(dot(caton,catocb)/(mag1*mag2));/* the N-Ca-Cb angle */

   switch(restype){
      case(ALA): /* just use the above ideal CB */
         break;
      case(ARG): /* arginine , assume trans (=PI) rotations for last two chi*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.ccc,geo.ccn,geo.cn,chiang[2],rotpos[NE]);
            cnext(geo,rotpos[CG],rotpos[CD],rotpos[NE],geo.ccn,geo.anc,geo.an,PI,rotpos[CZ]);
            cnext(geo,rotpos[CD],rotpos[NE],rotpos[CZ],geo.anc,geo.nan,geo.an,PI,rotpos[NH1]);
            mnext(rotpos[NE],rotpos[CZ],rotpos[NH1],geo.nan,geo.nan,geo.an,rotpos[NH2]);
         break;
      case(ASN):/* asparagine */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.acc,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.cao,geo.ao,chiang[1],rotpos[OD1]);
            mnext(rotpos[CB],rotpos[CG],rotpos[OD1],geo.cao,geo.can,geo.an,rotpos[ND2]);
         break;
      case(ASP):/* aspartate */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.acc,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.cao,geo.ao,chiang[1],rotpos[OD1]);
            mnext(rotpos[CB],rotpos[CG],rotpos[OD1],geo.cao,geo.cao,geo.ao,rotpos[OD2]);
         break;
      case(CYS):/* cysteine */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccs,geo.cs,chiang[0],rotpos[SG]);
         break;
      case(GLU):/* glutamate */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.acc,geo.ac,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.acc,geo.cao,geo.ao,chiang[2],rotpos[OE1]);
            mnext(rotpos[CG],rotpos[CD],rotpos[OE1],geo.cao,geo.cao,geo.ao,rotpos[OE2]);
         break;
      case(AA4):// artificial residue used to sample density
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.acc,geo.ac,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.acc,geo.cao,geo.ao,chiang[2],rotpos[OE1]);
//            mnext(rotpos[CG],rotpos[CD],rotpos[OE1],geo.cao,geo.cao,geo.ao,rotpos[OE2]);
         break;
      case(GLN):/* glutamine */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.acc,geo.ac,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.acc,geo.cao,geo.ao,chiang[2],rotpos[OE1]);
            mnext(rotpos[CG],rotpos[CD],rotpos[OE1],geo.cao,geo.can,geo.an,rotpos[NE2]);
         break;
      case(GLY):
         break;
      case(HIS):/*New - use variable ringang histadine  Note third chi = PI, rest of ring = 0*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,ringang,geo.bc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bcc,geo.cbn,geo.bn,chiang[1],rotpos[ND1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[ND1],geo.cbn,geo.bnb,geo.bn,PI,rotpos[CE1]);
            cnext(geo,rotpos[CG],rotpos[ND1],rotpos[CE1],geo.bnb,geo.nbn,geo.bn,0.0,rotpos[NE2]);
            cnext(geo,rotpos[ND1],rotpos[CE1],rotpos[NE2],geo.nbn,geo.bnb,geo.bn,0.0,rotpos[CD2]);
         break;
      case(ILE):// Major problem - tordyn assumes only terminal group has more than one atom
               //  Thus, treat all 4 terminal atoms as 1 fixed group - ASSUME chiang[i][1] = 169.0;
            chiang[1]= 169.0*PI/180.0;
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG1]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG1],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD1]);
#if 0
            mnext(rotpos[CA],rotpos[CB],rotpos[CG1],geo.ccc,ccc,geo.cc,rotpos[CG2]);
#endif
#if 1
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0]-tang,rotpos[CG2]);
#endif
         break;
      case(LEU):
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1]+tang,rotpos[CD2]);
         break;
      case(AA3): // artificial residue, used to sample density
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD1]);
//            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1]+tang,rotpos[CD2]);
         break;
      case(LYS): /* assume trans rot for last chi*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccc,geo.cc,chiang[1],rotpos[CD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD],geo.ccc,geo.ccc,geo.cc,chiang[2],rotpos[CE]);
            cnext(geo,rotpos[CG],rotpos[CD],rotpos[CE],geo.ccc,geo.ccn,geo.cn,PI,rotpos[NZ]);
         break;
      case(MET): /* assume trans rot for last chi*/
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.ccc,geo.ccs,geo.cs,chiang[1],rotpos[SD]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[SD],geo.ccs,geo.csc,geo.cs,chiang[2],rotpos[CE]);
         break;
      case(PHE): /*New - use variable ringang -  use aaa where want aac (not in table( */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,ringang,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.aaa,geo.aa,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD1],geo.aaa,geo.aaa,geo.aa,PI,rotpos[CE1]);
            cnext(geo,rotpos[CG],rotpos[CD1],rotpos[CE1],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CZ]);
            cnext(geo,rotpos[CD1],rotpos[CE1],rotpos[CZ],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CE2]);
            cnext(geo,rotpos[CE1],rotpos[CZ],rotpos[CE2],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CD2]);
         break; 
#if 0     
      case(PRO):
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.bbb,geo.bb,chiang[0],rotpos[CG]);
         cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bbb,geo.bbb,geo.bb,chiang[1],rotpos[CD]);
         break;
#endif
#if 1   // New = fix chiang  0   
      case(PRO):
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.bbb,geo.bb,0.0,rotpos[CG]);
         cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bbb,geo.bbb,geo.bb,0.0,rotpos[CD]);
         break;
#endif
      case(SER):
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccv,geo.cv,chiang[0],rotpos[OG]);
         break;
      case(THR):
           cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccv,geo.cv,chiang[0],rotpos[OG1]);
#if 0
           mnext(rotpos[CA],rotpos[CB],rotpos[OG1],geo.ccv,ccv,geo.cv,rotpos[CG2]);
#endif
           cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccv,geo.cv,chiang[0]+tang,rotpos[CG2]);
         break;
      case(TRP): /* New use variable ringang Note third chi = PI */
//            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.bcc,geo.bc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,ringang,geo.bc,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.bcc,geo.bbc,geo.bb,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD1],geo.bbc,geo.bbn,geo.bn,PI,rotpos[NE1]);
            cnext(geo,rotpos[CG],rotpos[CD1],rotpos[NE1],geo.bbn,geo.bng,geo.gn,0.0,rotpos[CE2]);
            cnext(geo,rotpos[CD1],rotpos[NE1],rotpos[CE2],geo.bng,geo.ggn,geo.gg,0.0,rotpos[CD2]);
            mnext(rotpos[CE2],rotpos[CD2],rotpos[CG],geo.bgg,geo.agg,geo.ag,rotpos[CE3]);
            cnext(geo,rotpos[CE2],rotpos[CD2],rotpos[CE3],geo.agg,geo.aag,geo.aa,0.0,rotpos[CZ3]);
            cnext(geo,rotpos[CD2],rotpos[CE3],rotpos[CZ3],geo.aag,geo.aaa,geo.aa,0.0,rotpos[CH2]);
            cnext(geo,rotpos[CE3],rotpos[CZ3],rotpos[CH2],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CZ2]);
         break;
      case(TYR): /* New use variable ringang  use aaa where want aac (not in table) */
            cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,ringang,geo.ac,chiang[0],rotpos[CG]);
            cnext(geo,rotpos[CA],rotpos[CB],rotpos[CG],geo.acc,geo.aaa,geo.aa,chiang[1],rotpos[CD1]);
            cnext(geo,rotpos[CB],rotpos[CG],rotpos[CD1],geo.aaa,geo.aaa,geo.aa,PI,rotpos[CE1]);
            cnext(geo,rotpos[CG],rotpos[CD1],rotpos[CE1],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CZ]);
            cnext(geo,rotpos[CD1],rotpos[CE1],rotpos[CZ],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CE2]);
            cnext(geo,rotpos[CE1],rotpos[CZ],rotpos[CE2],geo.aaa,geo.aaa,geo.aa,0.0,rotpos[CD2]);
            mnext(rotpos[CE1],rotpos[CZ],rotpos[CE2],geo.aaa,geo.aav,geo.av,rotpos[OH]);

         break;      
      case(VAL):
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0],rotpos[CG1]);
#if 0
         mnext(rotpos[CA],rotpos[CB],rotpos[CG1],geo.ccc,ccc,geo.cv,rotpos[CG2]);
#endif
         cnext(geo,rotpos[N],rotpos[CA],rotpos[CB],al1,geo.ccc,geo.cc,chiang[0]+tang,rotpos[CG2]);
         break;
            
    }
}


void setuprotomer(struct rotomer rot[])
// Reads rotlib, finds chiang and rotnum for each residue type and each rotomer
{
   int i;

   data_side(side,brdraw);
   for(i=ALA;i<=AA4;i++)
      rotlib(i,&rot[i].rotnum,rot[i].chiang); // returns library data for building all rotomers for ths rsesidue
}      


void checkrotset(struct griddata& den1,struct geometry& geo,int ifit,
        int resn,int *badset)
/*
   Checks for bad contacts between terminal atoms of residue
     (atomnum >=8) and set points
  geo.resside[NAFIT];// Top number of side chain atom - USED ONLY IF SIDE = RESIDUE
   checkatom returns 0 if point overlaps set
*/
{
   int i,geoloc;
    
    *badset = 0;
    geoloc = geo.nend[0]+resn;
    for(i=7;i<=geo.resside[geoloc];i++)
      if(checkatom(den1,geo,geo.set,geo.pos[geoloc][i])== 0){ // if atom bad
          *badset = *badset+1;
      }
}
        
void copygeotorotsav(struct geometry& geo,int resn,float rotsav[NAA][3])
{
   int i,k,geoloc;
 
   geoloc= geo.nend[0]+resn;
   for(i=5;i<=geo.resside[geoloc];i++){ // Copy each sidechain atom in geo
      for(k=0;k<=2;k++)
         rotsav[i][k] = geo.pos[geoloc][i][k];
   }
}
        
void copyrotsavtogeo(struct geometry& geo,int resn,float rotsav[NAA][3])
{
   int i,k,geoloc;
 
   geoloc= geo.nend[0]+resn;
   for(i=5;i<=geo.resside[geoloc];i++){ // Copy each sidechain atom in geo
      for(k=0;k<=2;k++)
        geo.pos[geoloc][i][k]= rotsav[i][k];
   }
}

/* Following stuff is required for Powell minimization:
   Powell routines are part of tracerecipes.inc which is called
     by trace.c++ -  Thus these routines should be considered as
     part of trace.  
   Also - the general Powell program is "powellmin", which is in trace.c++
     phead and pfunc are in trace.h which is included above.
   itranlim is a gloval variable, defined in tracae.c++ and used in powell routines
      to exit the routine
*/
extern int itranlim;  //Defined in trace.c++Used in tracerecipes = 0 if powell tran step exceeds allowed limit (must be called before include traacerecipes


void moveside(int iside,struct geometry& geo,struct tordyn& tordata,
          float ipos[NAFIT][NAA][3],float p[])
/*
Similar to onesidestep but used in Powelmin
  Adjustable parameters are pp[1,2,...nvar]   
    (substituted for  tordata.dqs[iunit-1]...     lunit in onsidestep
Copies initial pos (ipos) to tordata, rotates and
   returns new tordata
*/
{
    int i,k,sunit,lunit,count2,ii,jj,kk,iunit,i1,i2;
    float rvec[3],dist,rom[3][3],dchi,pt1[3],pt2[3];

    count2=0; /* counts the type 2 atoms, used to find N */
   i = iside;
    count2=count2+1;
    if(tordata.sidedata[i][0]!=0){
       sunit=tordata.sidedata[i][1]+1;/*start value ( sidedata[i][1]= Ca )*/
       lunit=tordata.sidedata[i][1]+tordata.sidedata[i][0]-1; /* last value of this sidechian */
       for(k=0;k<=2;k++){// Copy orig pos to tordata
          for(i1=sunit;i1<=lunit;i1++)
             tordata.uside[i1][0][k]=ipos[i1-sunit][0][k];
          i1=lunit;
          for(i2=1;i2<=tordata.sidedata[i][3];i2++)
            tordata.uside[i1][i2][k]=ipos[i1-sunit][i2][k];
       }//k
       for(k=0;k<=2;k++)  /* copy Cb to uside */
               tordata.uside[sunit][0][k]=tordata.unit[0][i][1][k];
       if(geo.dir==0) // N to C direciton
              ii = count2; 
       else if(geo.dir==1)
              ii = count2+1;
       for(k=0;k<=2;k++) //first rotation vector = Ca to Cb
             pt1[k] = tordata.unit[0][i][0][k];  // CA
       for(ii= 1;ii< tordata.sidedata[i][0];ii++){
            iunit = sunit+ii -1;
            for(k=0;k<=2;k++){ //first rotation vector = Ca to Cb
               pt2[k] = tordata.uside[iunit][0][k];
               rvec[k] = pt2[k]-pt1[k];
            }
            dist = sqrt( dot(rvec,rvec) );
            for(k=0;k<=2;k++) 
               rvec[k] = rvec[k]/dist;
            dchi = p[ii];
            rotaxis(geo.de,geo.ep,dchi,rvec,rom); /* rom is rotation matrix for rotation by dchi about unit vector rvec which is assumed to be at origin*/  
            for(kk=iunit+1;kk<=lunit;kk++) // Rotate all subsequent single pts
               rotateset(pt2,tordata.uside[kk][0],rom); // pt2 is origin
            for(jj=1;jj<=tordata.sidedata[i][3];jj++) // rotate all points in last unit connected to connect point
                  rotateset(pt2,tordata.uside[lunit][jj],rom); // pt2 is origin
            for(k=0;k<=2;k++)
               pt1[k] = pt2[k]; 
      }//ii
   } /* if sidedata */
}                                


void movesidecb(int iside,struct geometry& geo,struct tordyn& tordata,
          float ipos[NAFIT][NAA][3],float p[])
/*
Modidy moveside so that allows rigid translation of cb (or Ca and rest of side chain.
  Cb[x,y,z] translation are in p[chimax,chimax+1,chimax+2]/4
     where chimax = tordata.sidedata[i][0];
  NOTE:  divide by 4 because initial step size set = 1
Similar to onesidestep but used in Powelmin
  Adjustable parameters are pp[0,1,2,...nvar]   
    (substituted for  tordata.dqs[iunit-1]...     lunit in onsidestep
Copies initial pos (ipos) to tordata, rotates and
   returns new tordata
Two possibilities:
  1) - Translate both Ca and Cb (and all other side atoms)  by same amount so direction of Ca-Cb vector does not change
  2) - Translate just Cb (and all other side atoms) -so that Ca=Cb vector changes direction
Doesn't seem to make much difference.
NOTE: THIS PROCESS IS NON-LINEar AND DOES NOT CONVERGE TO  BEST ANSWER -
   But gives better answer each time called.
  because, initially, translate all side chains by delta Cb.
Next time called, starting with whole side chain rotated so
   corresponds to new Ca-Cb bond, so get different result.
If Call 4 or 5 times in a row - roughly converges to best value.  
*/
{
    int i,k,sunit,lunit,count2,ii,jj,kk,iunit,i1,i2;
    float rvec[3],dist,rom[3][3],dchi,pt1[3],pt2[3];
    int chimax;

    count2=0; /* counts the type 2 atoms, used to find N */
   i = iside;
    count2=count2+1;
    if(tordata.sidedata[i][0]!=0){
       chimax = tordata.sidedata[i][0];
       sunit=tordata.sidedata[i][1]+1;/*start value ( sidedata[i][1]= Ca )*/
       lunit=tordata.sidedata[i][1]+tordata.sidedata[i][0]-1; /* last value of this sidechian */
       for(k=0;k<=2;k++){// Copy orig pos to tordata and add Cb translation
          for(i1=sunit;i1<=lunit;i1++)
             tordata.uside[i1][0][k]=ipos[i1-sunit][0][k]+p[chimax+k]/4.0;
          i1=lunit;
          for(i2=1;i2<=tordata.sidedata[i][3];i2++)
            tordata.uside[i1][i2][k]=ipos[i1-sunit][i2][k]+p[chimax+k]/4.0;
       }//k
       for(k=0;k<=2;k++)  /* copy Cb to uside and add varible step*/
               tordata.uside[sunit][0][k]=tordata.unit[0][i][1][k]+p[chimax+k]/4.0;
       if(geo.dir==0) // N to C direciton
              ii = count2; 
       else if(geo.dir==1)
              ii = count2+1;
#if 0 //new 14dec00 = try  translate Ca same as Cb Does not seem as good - does not change plane of rings?? 
       for(k=0;k<=2;k++) //first rotation vector = Ca to Cb
             pt1[k] = tordata.unit[0][i][0][k]+p[chimax+k]/4.0;  // CA
#endif
#if 1 // Do not translate Ca 
       for(k=0;k<=2;k++) //first rotation vector = Ca to Cb
             pt1[k] = tordata.unit[0][i][0][k];  // CA
#endif
       for(ii= 1;ii< tordata.sidedata[i][0];ii++){
            iunit = sunit+ii -1;
            for(k=0;k<=2;k++){ //first rotation vector = Ca to Cb
               pt2[k] = tordata.uside[iunit][0][k];
               rvec[k] = pt2[k]-pt1[k];
            }
            dist = sqrt( dot(rvec,rvec) );
            for(k=0;k<=2;k++) 
               rvec[k] = rvec[k]/dist;
            dchi = p[ii];
            rotaxis(geo.de,geo.ep,dchi,rvec,rom); /* rom is rotation matrix for rotation by dchi about unit vector rvec which is assumed to be at origin*/  
            for(kk=iunit+1;kk<=lunit;kk++) // Rotate all subsequent single pts
               rotateset(pt2,tordata.uside[kk][0],rom); // pt2 is origin
            for(jj=1;jj<=tordata.sidedata[i][3];jj++) // rotate all points in last unit connected to connect point
                  rotateset(pt1,tordata.uside[lunit][jj],rom); // pt1 is origin
            for(k=0;k<=2;k++)
               pt1[k] = pt2[k]; 
      }//ii
   } /* if sidedata */
}                                

void movesidecbca(int iside,struct geometry& geo,struct tordyn& tordata,
          float ipos[NAFIT][NAA][3],float p[])
/*
Modidy moveside so that allows rigid translation of cb (or Ca and rest of side chain.
  Cb[x,y,z] translation are in p[chimax,chimax+1,chimax+2]/4
     where chimax = tordata.sidedata[i][0];
  NOTE:  divide by 4 because initial step size set = 1
Similar to onesidestep but used in Powelmin
  Adjustable parameters are pp[0,1,2,...nvar]   
    (substituted for  tordata.dqs[iunit-1]...     lunit in onsidestep
Copies initial pos (ipos) to tordata, rotates and
   returns new tordata
Two possibilities:
  1) - Translate both Ca and Cb (and all other side atoms)  by same amount so direction of Ca-Cb vector does not change
  2) - Translate just Cb (and all other side atoms) -so that Ca=Cb vector changes direction
Doesn't seem to make much difference.
NOTE: THIS PROCESS IS NON-LINEar AND DOES NOT CONVERGE TO  BEST ANSWER -
   But gives better answer each time called.
  because, initially, translate all side chains by delta Cb.
Next time called, starting with whole side chain rotated so
   corresponds to new Ca-Cb bond, so get different result.
If Call 4 or 5 times in a row - roughly converges to best value. 
New: Crude way to allow Ca-Cb vector to have arbitrary orientation:
  Cb[x,y,z] translation are in p[chimax,chimax+1,chimax+2]/4
  Ca[x,y,z] translation are in p[chimax+3,chimax+4,chimax+5]/4
  1) Save initial Ca-Cb vector = rvec0 
  2) Translate Cb by variable amount Cb'.
  4) Find new Ca-Cb' vector = rvec.
  5) Make vector rrot = rvec0 X rvec which is Perpendicualr to rvec0 and rvec.
  6) Rotatate initial Cb and rest of side about rrot by 
         angle phi = angle  between rvec0 and rvec abount orig Ca.
  7) Then translate all atoms by Ca tran
*/
{
    int i,k,sunit,lunit,count2,ii,jj,kk,iunit,i1,i2;
    float rvec0[3],rrot[3],rvec[3],dist,rom[3][3],dchi,pt1[3],pt2[3];
    float ang,cavec[3],cbvec[3],cb2vec[3];
    int chimax;

    count2=0; /* counts the type 2 atoms, used to find N */
    i = iside;
    count2=count2+1;
    if(tordata.sidedata[i][0]!=0){
       chimax = tordata.sidedata[i][0];
       sunit=tordata.sidedata[i][1]+1;/*start value ( sidedata[i][1]= Ca )*/
       lunit=tordata.sidedata[i][1]+tordata.sidedata[i][0]-1; /* last value of this sidechian */
       for(k=0;k<=2;k++){ // initial Ca  and Cbvector
          cavec[k]=tordata.unit[0][i][0][k];
          cbvec[k]=tordata.unit[0][i][1][k];
       }
       for(k=0;k<=2;k++) // initial vector = Ca to Cb
          rvec0[k]= cbvec[k]- cavec[k];
       dist = sqrt( dot(rvec0,rvec0) );
       for(k=0;k<=2;k++) 
           rvec0[k] = rvec0[k]/dist;
       for(k=0;k<=2;k++) // translate Cb by cb tran
          cb2vec[k] = cbvec[k]+p[chimax+k]/4.0;
       for(k=0;k<=2;k++) // new vector = Ca to Cb
          rvec[k]= cb2vec[k]-cavec[k];
       dist = sqrt( dot(rvec,rvec) );
       for(k=0;k<=2;k++) 
               rvec[k] = rvec[k]/dist;
       crossprod(rvec0,rvec,rrot); // rrot is perp. to rvec0, and rvec
       dist = sqrt( dot(rrot,rrot) );
       for(k=0;k<=2;k++) 
           rrot[k] = rrot[k]/dist;
       ang = acos( dot(rvec0,rvec) ); // angle between rev0, and rvec
       rotaxis(geo.de,geo.ep,ang,rrot,rom); /* rom is rotation matrix for rotation by dchi about unit vector rvec which is assumed to be at origin*/  

       for(k=0;k<=2;k++){// Copy orig pos to tordata
          for(i1=sunit;i1<=lunit;i1++)
             tordata.uside[i1][0][k]=ipos[i1-sunit][0][k];
          i1=lunit;
          for(i2=1;i2<=tordata.sidedata[i][3];i2++)
            tordata.uside[i1][i2][k]=ipos[i1-sunit][i2][k];
       }//k

       for(kk=sunit;kk<=lunit;kk++) // Rotate all subsequent single pts about origin = cavec
              rotateset(cavec,tordata.uside[kk][0],rom); //cavec is pt on rotation axis is origin
       for(jj=1;jj<=tordata.sidedata[i][3];jj++) // rotate all points in last unit connected to connect point
              rotateset(cavec,tordata.uside[lunit][jj],rom); // pt1 is origin
       
       for(k=0;k<=2;k++){// add Ca translation to all pts
          for(i1=sunit;i1<=lunit;i1++)
             tordata.uside[i1][0][k]=tordata.uside[i1][0][k]+p[chimax+k+3]/4.0;
          i1=lunit;
          for(i2=1;i2<=tordata.sidedata[i][3];i2++)
            tordata.uside[i1][i2][k]=tordata.uside[i1][i2][k]+p[chimax+k+3]/4.0;
       }//k
       for(k=0;k<=2;k++) // add Ca translation to Ca to make new Ca = pt1
             pt1[k] = cavec[k]+p[chimax+k+3]/4.0;  // CA

       for(ii= 1;ii< tordata.sidedata[i][0];ii++){ // Now rotate about bonds
            iunit = sunit+ii -1;
            for(k=0;k<=2;k++){ //first rotation vector = Ca to Cb
               pt2[k] = tordata.uside[iunit][0][k];
               rvec[k] = pt2[k]-pt1[k];
            }
            dist = sqrt( dot(rvec,rvec) );
            for(k=0;k<=2;k++) 
               rvec[k] = rvec[k]/dist;
            dchi = p[ii];
            rotaxis(geo.de,geo.ep,dchi,rvec,rom); /* rom is rotation matrix for rotation by dchi about unit vector rvec which is assumed to be at origin*/  
            for(kk=iunit+1;kk<=lunit;kk++) // Rotate all subsequent single pts
               rotateset(pt2,tordata.uside[kk][0],rom); // pt2 is origin
            for(jj=1;jj<=tordata.sidedata[i][3];jj++) // rotate all points in last unit connected to connect point
                  rotateset(pt1,tordata.uside[lunit][jj],rom); // pt1 is origin
            for(k=0;k<=2;k++)
               pt1[k] = pt2[k]; 
      }//ii
   } /* if sidedata */
}                                

float rotfuncdencb(float p[],phead)
/* 
New 9feb01 - determines mindist of far side chain atoms from mainchain -
   If < geo.sidemaindist - then increase nbad
Modifiy rotfuncden so that allows Cb to move.
  Cb[x,y,z] translation are in p[chimax+1,chimax+2,chimax=3]
     where chimax = tordata.sidedata[i][0];
  Check for total distance Cb moved, if greater than (say) 0.5 A,
    , then increase nbad by 1, this should limit movment.
   
Returns the side residue quality func
*/
{
    int nbad,badset;
    float aveden,value,cbdist;
    float mindist;

//    maxdist = 0.5; //Now replaced by geo.maxcbmove   maximum allowed movementof Cb
    itranlim=1;  // dummy variable  - if = 0, then powel quits;
    movesidecb(1,geo,tordata,ipos,p); //iside = 1
    avedensetoneside(den1,geo,tordata,1,&aveden,&nbad,&badset);
//    avestddevoneside(den1,geo,tordata,iside,&aveden,&*nbad,&badset,&stddev); 
    if(geo.quit==1){ // means atom out of range
         return 10000.0;  // very bad value
    }
#if 1 // New 9feb01 - find mindist of far atoms from mainchain geo atoms - and modify badset
//geo.sidemaindist=2.5;
    mindist = 1000.0;
    findside_maindisttor(geo,geo.resn,1,&mindist);
    if(mindist < 900.0){ // Means there were atoms to be tested
       if(mindist< geo.sidemaindist)
          badset = badset+1;
    }
#endif
    cbdist =  sqr(p[nvar]) + sqr(p[nvar-1])+sqr(p[nvar-2])/16.0; //divide by 16 because step size = p/4 = dist from initial Cb moved
    if(cbdist > sqr(geo.maxcbmove) ){ // Max allowed movement of 0.5? angstroms
       nbad = nbad+3;
//printf("cbdist = %5.3f\n",sqrt(cbdist) );
    }
    value = -(aveden - maidfile.denscale*(geo.valueweight[0]*nbad + geo.valueweight[1]*badset) );
//    value = -(aveden - maidfile.denscale*(0.5*nbad + 5.0*badset) );  //change sign becuse want maximum value, but Powell finds minimum 
//printf(" value = %5.3f p[1] = %5.3f  p[2]= %5.3f\n",value,p[1],p[2]);
    return value;
}   

float rotfuncdencbnoset(float p[],phead)
/* 
Modifiy rotfuncden so that allows Cb to move.
  Cb[x,y,z] translation are in p[chimax+1,chimax+2,chimax=3]
     where chimax = tordata.sidedata[i][0];
  Check for total distance Cb moved, if greater than (say) 0.5 A,
    , then increase nbad by 1, this should limit movment.
   
Returns the side residue quality func
*/
{
    int nbad;
    float aveden,value,cbdist;

//    maxdist = 0.5; //Now replaced by geo.maxcbmove   maximum allowed movementof Cb
    itranlim=1;  // dummy variable  - if = 0, then powel quits;
    movesidecb(1,geo,tordata,ipos,p); //iside = 1
    avedenoneside(den1,geo,tordata,1,&aveden,&nbad); // find density for side chain = 1
    cbdist =  sqr(p[nvar]) + sqr(p[nvar-1])+sqr(p[nvar-2])/16.0; //divide by 16 because step size = p/4 = dist from initial Cb moved
    if(cbdist > sqr(geo.maxcbmove) ){ // Max allowed movement of 0.5 angstroms
        nbad = nbad +5;
    }
    value = -(aveden - maidfile.denscale*(geo.valueweight[0]*nbad));
//    value = -(aveden - maidfile.denscale*0.5*nbad );  //change sign becuse want maximum value, but Powell finds minimum 
//printf(" value = %5.3f p[1] = %5.3f  p[2]= %5.3f\n",value,p[1],p[2]);
    return value;
}   

float rotfuncden(float p[],phead)
/* Returns the side residue quality func
*/
{
    int nbad,badset;
    float aveden,value;

    itranlim=1;  // dummy variable  - if = 0, then powel quits;
#if 0
for(i=1;i<=nvar;i++)
 p[i]=p[i]/10;
#endif
    moveside(1,geo,tordata,ipos,p); //iside = 1
    avedensetoneside(den1,geo,tordata,1,&aveden,&nbad,&badset); 
    if(geo.quit==1){ // means atom out of range
         return 10000.0;  // very bad value
    }
    value = -(aveden - maidfile.denscale*(geo.valueweight[0]*nbad + geo.valueweight[1]*badset) );
//    value = -(aveden - maidfile.denscale*(0.5*nbad + 5.0*badset) );  //change sign becuse want maximum value, but Powell finds minimum 
#if 0
printf("value = %5.3f",value);
for(i =1;i<=nvar;i++)
   printf("  p[%d] = %5.3f",i,p[i]);
printf("\n");
#endif
    return value;
}
   
float rotfuncdennoset(float p[],phead)
/* Returns the side residue quality func
*/
{
    int nbad;
    float aveden,value;

    itranlim=1;  // dummy variable  - if = 0, then powel quits;
    moveside(1,geo,tordata,ipos,p); //iside = 1
    avedenoneside(den1,geo,tordata,1,&aveden,&nbad); // find density for side chain = 1
    value = -(aveden - maidfile.denscale*(geo.valueweight[0]*nbad) );
//    value = -(aveden - maidfile.denscale*0.5*nbad );  //change sign becuse want maximum value, but Powell finds minimum 
//printf(" value = %5.3f p[1] = %5.3f  p[2]= %5.3f\n",value,p[1],p[2]);
    return value;
}   

void bestrotpow(int useset,int movecb,struct griddata& den1,struct geometry& geo,float *aveden,
        int *nbad,float pp[])
/* Use Powell minimiazation to maximize density
NOTE:  nbad is no longer found - returned value is nonsense.
  If movecb = 1, allow cb to move,if = 0 , keep cb fixed
  If movecb = 2 allow both Ca and Cb to move // NOTE: this did not work
     is now deleted - see assignseq.c27oct99
*/
{
     int i,i1,k,i2,nvar,lunit,sunit,iside;
// Following are not used but are passed in pfunc
     struct pepfit fit1,fit2; //These are not used but are needed in pfunc
     float origphi[3],origpsi[3],ipos[NAFIT][NAA][3],fpos[NAFIT][NAA][3];
     float zaxis[3],xaxis[3]; // this is needed for pfunc definition but is not used.

     int cafix =0; // this is needed for pfunc definition but is not used.
      iside = geo.oneside-1;
      geo.dir = 0;// N to C direction
      geo.begend[0]= geo.nend[0]+iside;
      geo.begend[1]=geo.begend[0]+2;
     ncopygeotoauto(geo,tordata); 
#if 0 // debug - check avedengeooneside - OK BUT Counts some atoms twice
float oneden;
int onebad;
    avedenoneside(den1,geo,tordata,1,&oneden,&onebad); //NEED to find nbad 
    printf(" for tordata: oneden = %5.3f, onebad = %d\n",oneden,onebad);
    avedengeooneside(den1,geo,geo.pos,geo.nend[0]+1,&oneden,&onebad);
    printf(" for geo: oneden = %5.3f, onebad = %d\n",oneden,onebad);
    exit(0);
#endif
     i = 1;  // = iside
     sunit=tordata.sidedata[i][1]+1;/*start value ( sidedata[i][1]= Ca )*/
     lunit=tordata.sidedata[i][1]+tordata.sidedata[i][0]-1; /* last value of this sidechian */
     for(k=0;k<=2;k++){// Store original positions in ipos
          for(i1=sunit;i1<=lunit;i1++)
             ipos[i1-sunit][0][k]=tordata.uside[i1][0][k];
          i1=lunit;
          for(i2=1;i2<=tordata.sidedata[i][3];i2++)
            ipos[i1-sunit][i2][k]=tordata.uside[i1][i2][k];
     }//k
     if(movecb ==0){ // fix cb
        nvar = lunit-sunit+1; //=(tordata.sidedata[i][0]-1) variable run from 1, .. nvar
        if(useset ==1)
           powellmin(aveden,pfunc,rotfuncden); // aveden is really a quality func - see value in rotfuncden
        else
           powellmin(aveden,pfunc,rotfuncdennoset); // aveden is really a quality func - see value in rotfuncden
        moveside(1,geo,tordata,ipos,pp); //iside = 1
     }
     else if(movecb ==1){
           nvar = lunit-sunit+1+3; //=(tordata.sidedata[i][0]-1) variable run from 1, .. nvar
        if(useset ==1)
           powellmin(aveden,pfunc,rotfuncdencb); // aveden is really a quality func - see value in rotfuncden
        else
           powellmin(aveden,pfunc,rotfuncdencbnoset); // aveden is really a quality func - see value in rotfuncden
        movesidecb(1,geo,tordata,ipos,pp); //iside = 1
#if 0 
float cbdist = sqrt( sqr(pp[nvar]) + sqr(pp[nvar-1])+sqr(pp[nvar-2]) )/4.0; //divide by 4 because step size = p/4  = dist from initial Cb moved
printf(" Cbdist = %5.3f\n",cbdist);
#endif
     }
//printf("At end aveden = %f  fden = %f\n",*aveden,fden);
}

void refineonerotgeo(int useset,struct griddata& den1,struct geometry& geo,
          struct rotomer rot[],int resn,int rottype,
          float *bestqual,int *bestnbad,int *bestbadset,int *bestrotn)
/*
NOTE: THIS USES geo.movecb in place of movecb
Indentical to refineonerot except uses
  buildrot(geo,geo.pos[resn][4],geo.pos[resn][1],
                  geo.pos[resn][5],rottype,rot[rottype].chiang[rotn],rotpos);
  (In refineonerotgeo resn is the the geo number, while in refineonerot()
         it is the fit number)
    set geo.oneside = resn - geo.nend[0]
    and use resn in place of geo.onesav
    and resn-geo.nend[0] in place of geo.onesav
    and sets geoloc = resn.
*/
{
    int quitcond,rotn,nbad,badset,usepowell;
    int rotbeg,rotend;
    float rotpos[NATOMS][3];
    float rotsav[NAA][3],qual,bestden,pp[NDIM];
    int origoneside;

#if 0 // New - 19jun00 - now defined globally in geo.valueweight
    valueweight[0] = 0.5; // value = -(aveden - maidfile.denscale*(valueweight[0]*nbad + valueweight[1]*badset) );
    valueweight[1] = 1.0;
#endif
    origoneside = geo.oneside;
    geo.oneside = resn- geo.nend[0];
    usepowell = 1; //NOTE: do not set = 0, NOT DEBUGGED!if = -1, do not refine, if = 0 use tordynm,  if = 1 then use powell min,else
    if(usepowell == -1)
      printf(" NO REFINEMENT\n");
    *bestqual = -1000.0;
    if( *bestrotn != -1){
        rotbeg = rotend = *bestrotn;
    }
    else{ // try all rotomers
       rotbeg = 1;
       rotend = rot[rottype].rotnum;
    }
    for(rotn=rotbeg;rotn<=rotend;rotn++){// try each rotomer for rottype 
//printf("rottype = %d, rotn = %d\n",rottype,rotn);
      geo.resn=geo.begend[0]+1; //Only copy 1res on either side to geo - this is number of geo that side is attached to
      buildrot(geo,geo.pos[resn][4],geo.pos[resn][1],
                  geo.pos[resn][5],rottype,rot[rottype].chiang[rotn],rotpos);
      copyonerottogeo(rotpos,geo,resn,rottype);// copy rotomer side chain to geo
      if((rottype==GLY)||(rottype==ALA)||(rottype==PRO) ){// do not need dynamics
//      if((rottype==GLY)||(rottype==ALA)||(rottype==GLY ) ){
          geo.dir = 0;// N to C direction
          geo.begend[0]= geo.nend[0]+geo.oneside-1;
          geo.begend[1]=geo.begend[0]+2;
          ncopygeotoauto(geo,tordata); 
//          avedenoneside(den1,geo,tordata,1,&bestden,&nbad); //NEED to find nbad and true aveden
          avedensetoneside(den1,geo,tordata,1,&bestden,&nbad,&badset); 
          *bestrotn = rotn;
          *bestnbad = nbad;
          *bestbadset = badset;
          *bestqual = bestden - maidfile.denscale*geo.valueweight[0]*nbad ;  //Do not consider badset for these(alway bad for PRO)
          goto skipdyn;
      }
      if(usepowell==0) // This is not debugged (NOTE> quicond not defined
          tordynbat(quitcond,den1,geo,tordata,&bestden,&nbad);// Run torsional dynamics on side chain, geo.fixend is not used
      else if(usepowell == 1) { //Do not refine proline;  Use Powell minimization to find best fit
          bestrotpow(useset,geo.movecb,den1,geo,&bestden,&nbad,pp); // nbad is nonsense -bestden is a qual function - not a pure density
          copyonesidetogeo(geo,tordata);
       }
//bestden = 50.0;
      qual = bestden;
      if(qual> *bestqual){
          avedensetoneside(den1,geo,tordata,1,&bestden,&nbad,&badset); 
          *bestqual= qual;
          *bestbadset = badset;
          *bestrotn = rotn;
          *bestnbad = nbad;
          copygeotorotsav(geo,resn-geo.nend[0],rotsav);
       }
    } //rotn;
#if 0
printf("For ifit = %d resn= %d  rottype = %d   Best rot = %d: bestnbad = %d bestqual = %5.3f\n",
        ifit,resn,rottype,*bestrotn,*bestnbad,*bestqual);
#endif
    copyrotsavtogeo(geo,resn-geo.nend[0],rotsav); 
  skipdyn:; // jump to here if only 1 rotomer and no dynamics
  geo.oneside= origoneside;
}

void refineonerot(int useset,struct griddata& den1,struct geometry& geo,
          struct rotomer rot[],int ifit,int resn,int rottype,
          float *bestqual,int *bestnbad,int *bestbadset,int *bestrotn)
/*
New: if bestrotn = -1; then try all rotomers, else just use bestrotn
      Used when call routine several times in row.
 Builds rotomer, copies to geo and refines position using dynamcis
  Returns the best value and best rotomer as a geo file
    quality defined  by:
      Old def.   qual = bestden * (1.0 - 0.3 nbad - 0.3 badset); (No good because bestden may be < 0
      New def.   qual  = bestden - 0.3*nbad*maidfile.denscale - 0.3*badset*maidfile.denscale
*/
{
    int quitcond,rotn,geoloc,nbad,badset,usepowell;
    int rotbeg,rotend;
    float rotpos[NATOMS][3];
    float rotsav[NAA][3],qual,bestden,pp[NDIM];
    float pheorig,tyrorig;

    usepowell = 1; //NOTE: do not set = 0, NOT DEBUGGED!if = -1, do not refine, if = 0 use tordynm,  if = 1 then use powell min,else
    if(usepowell == -1)
      printf(" NO REFINEMENT\n");
    pheorig = geo.phe; // angle of ring in phenylalanine
    tyrorig = geo.tyr;
    geoloc=geo.nend[0]+geo.oneside;
    *bestqual = -1000.0;
    if( *bestrotn != -1){
        rotbeg = rotend = *bestrotn;
    }
    else{ // try all rotomers
       rotbeg = 1;
       rotend = rot[rottype].rotnum;
    }
//rotend = 5;rotbeg = 5;
  float delphe,phe[4],tyr[4];
  int iphe,nphe;
  delphe = 7.0*PI/180.0;  // delta phenylalanine angle = 3 degrees
  nphe = 1;  // default
  if(geo.ringang==1){
   if( (rottype == PHE)||(rottype==HIS)||(rottype==TRP) ){
      nphe = 3;
      phe[1]=pheorig;
      phe[2] = pheorig-delphe;
      phe[3] = pheorig+delphe;
   }
   else if(rottype==TYR){ 
      nphe = 3;
      tyr[1]=tyrorig;
      tyr[2] = tyrorig-delphe;
      tyr[3] = tyrorig+delphe;
   }
  }
  for(iphe=1;iphe<=nphe;iphe++){ // if ring, try different geo.phe angles
    if(geo.ringang==1){
       if( (rottype == PHE)||(rottype==HIS)||(rottype==TRP) )
          geo.phe = phe[iphe];
       else if(rottype==TYR)
          geo.tyr = tyr[iphe];
    }
    for(rotn=rotbeg;rotn<=rotend;rotn++){// try each rotomer for rottype 
      if(fit[ifit].nside[resn] == -GLY) // must add Cb
          findcb(geo,fit[ifit].pos[resn][4],fit[ifit].pos[resn][1],
               fit[ifit].pos[resn][3],fit[ifit].pos[resn][5]);          
      geo.resn=geo.begend[0]+1; //(See copyonerottogeo()  - Only copy 1res on either side to geo - this is number of geo that side is attached to
      buildrot(geo,fit[ifit].pos[resn][4],fit[ifit].pos[resn][1],
                  fit[ifit].pos[resn][5],rottype,rot[rottype].chiang[rotn],rotpos);
      copyonerottogeo(rotpos,geo,geoloc,rottype);// copy rotomer side chain to geo
      if((rottype==GLY)||(rottype==ALA)||(rottype==PRO) ){// do not need dynamics
//      if((rottype==GLY)||(rottype==ALA)||(rottype==GLY ) ){
          geo.dir = 0;// N to C direction
          geo.begend[0]= geo.nend[0]+geo.oneside-1;
          geo.begend[1]=geo.begend[0]+2;
          ncopygeotoauto(geo,tordata); 
//          avedenoneside(den1,geo,tordata,1,&bestden,&nbad); //NEED to find nbad and true aveden
          avedensetoneside(den1,geo,tordata,1,&bestden,&nbad,&badset); 
          *bestnbad = nbad;
          *bestbadset = badset;
          *bestqual = bestden - maidfile.denscale*geo.valueweight[0]*nbad ;  //Do not consider badset for these(alway bad for PRO)
          goto skipdyn;
      }
      if(usepowell==0)
          tordynbat(quitcond,den1,geo,tordata,&bestden,&nbad);// Run torsional dynamics on side chain, geo.fixend is not used
      else if(usepowell == 1) { //Do not refine proline;  Use Powell minimization to find best fit
#if 0 // New 18may00 - use geo.movecb
          movecb = 0; // keep Cb position fixed
//          movecb = 1; // move Cb
#endif
//printf("geo.movecb = %d\n",geo.movecb);
          bestrotpow(useset,geo.movecb,den1,geo,&bestden,&nbad,pp); // bestden is a qual function - not a pure density
          if(bestden< -900.0){
            if(verbose>=1)
               printf("ERROR: fit %d  res %d rotn %d: atom out of density in bestrotpow - using intial default rotomer\n",ifit,resn,rotn);
            *bestqual = -1000.0;
            goto skipdyn;
          }
          else
             copyonesidetogeo(geo,tordata);
     }
//      checkrotset(den1,geo,ifit,geo.oneside,&badset); // Not needed because checked in bestrogpow
//printf("badset = %d\n",badset);
//       qual = (*bestden)*(1.0 - 0.3*(*nbad) - 0.3*(*badset) ); // Old def
//       qual = (bestden) - 0.3*(nbad)*maidfile.denscale - 0.3*(badset)*maidfile.denscale;
//       qual = (bestden) - 0.5*(nbad)*maidfile.denscale - 5.0*(badset)*maidfile.denscale;
         qual = bestden;
#if 0
if(rottype==ARG)   
  printf(" For ifit =%d resn =%d rottype=%d  rotn = %d:   bestden = %5.3f  nbad = %d  badset = %d\n",
     ifit,resn,rottype,rotn,bestden,nbad,badset);    
#endif
      if(qual> *bestqual){
//          *bestnbad = nbad;// nbad is no longer found in bestrotpow
          *bestqual= qual;
          *bestbadset = badset;
          *bestrotn = rotn;
          if(usepowell==0)
             copygeotorotsav(geo,resn,rotsav);
          else if(usepowell == 1){ // last tordata residue fit
            copygeotorotsav(geo,geo.oneside,rotsav); //geo.oneside = resn
          }
       }
//printf("qual = %5.3f  rottype = %d, rotn = %d\n",qual,rottype,rotn);
    } //rotn;
  }// iphe
#if 0 // Use this rotomer and allow Cb to move
    copyrotsavtogeo(geo,geo.oneside,rotsav); 
   geo.movecb = 1; // allow Cb to move
    bestrotpow(useset,geo.movecb,den1,geo,&bestden,&nbad,pp); // bestden is a qual function - not a pure density
    *bestqual = bestden;
    copyonesidetogeo(geo,tordata);
    copygeotorotsav(geo,geo.oneside,rotsav); //geo.oneside = resn
    geo.movecb = 0;
#endif
#if 0
printf("For ifit = %d resn= %d  rottype = %d   Best rot = %d: bestnbad = %d bestqual = %5.3f\n",
        ifit,resn,rottype,*bestrotn,*bestnbad,*bestqual);
#endif
    copyrotsavtogeo(geo,geo.oneside,rotsav); 
  skipdyn:; // jump to here if only 1 rotomer and no dynamics
#if 1 // This is needed to to return correct bestnbad
 ncopygeotoauto(geo,tordata); 
 avedenoneside(den1,geo,tordata,1,&bestden,bestnbad); // find true density for side chain = 1 - 
//printf("bestden = %5.3f  bestnbad = %d\n",bestden,*bestnbad); // bestden same as above
#endif
 geo.phe = pheorig;
 geo.tyr=tyrorig;
}


void addoneletter(char nm[1],int *resn)
{
  if(strcasecmp(nm,"A")==0) *resn=ALA;
     else if(strcasecmp(nm,"R")==0) *resn=ARG;
     else if(strcasecmp(nm,"N")==0) *resn=ASN;
     else if(strcasecmp(nm,"D")==0) *resn=ASP;
     else if(strcasecmp(nm,"C")==0) *resn=CYS;
     else if(strcasecmp(nm,"Q")==0) *resn=GLN;
     else if(strcasecmp(nm,"E")==0) *resn=GLU;
     else if(strcasecmp(nm,"G")==0) *resn=GLY;
     else if(strcasecmp(nm,"H")==0) *resn=HIS;
     else if(strcasecmp(nm,"I")==0) *resn=ILE;
     else if(strcasecmp(nm,"L")==0) *resn=LEU;
     else if(strcasecmp(nm,"K")==0) *resn=LYS;
     else if(strcasecmp(nm,"M")==0) *resn=MET;
     else if(strcasecmp(nm,"F")==0) *resn=PHE;
     else if(strcasecmp(nm,"P")==0) *resn=PRO;
     else if(strcasecmp(nm,"S")==0) *resn=SER;
     else if(strcasecmp(nm,"T")==0) *resn=THR;
     else if(strcasecmp(nm,"W")==0) *resn=TRP;
     else if(strcasecmp(nm,"Y")==0) *resn=TYR;
     else if(strcasecmp(nm,"V")==0) *resn=VAL;
     else {
      printf("ERROR:  %s is not in  residue dictionary\n",nm);
      exit(0);
     }
}

void countside(struct geometry& geo,struct tordyn& tordata,
      struct rotomer rot[],int sidecount[]) 
/* similar to avedenoneside but just counts atoms
   in side chain after convert to  tordata format
*/
{
   int i,j,ns,count,iside,ifit,resn,rottype,rotn,geoloc;
   float rotpos[NATOMS][3];

    iside = 1;
    ifit = 1;
    resn=1; 
    rotn = 1;
    copyregfittogeo(geo,fit[ifit],resn-1,resn+1);
    geo.dir = 0;// N to C direction
    geo.begend[0]= geo.nend[0]+resn-1;
    geo.begend[1]=geo.begend[0]+2;
    geoloc = geo.nend[0] + resn;
    for(rottype = ALA; rottype<= AA4; rottype++){
      buildrot(geo,fit[ifit].pos[resn][4],fit[ifit].pos[resn][1],
          fit[ifit].pos[resn][5],rottype,rot[rottype].chiang[rotn],rotpos);
      copyonerottogeo(rotpos,geo,geoloc,rottype);// copy rotomer side chain to geo
      ncopygeotoauto(geo,tordata); 
      ns=tordata.sidedata[iside][1];// start number for this chain = Ca
      count = 0;
      if(tordata.sidedata[iside][0]==0) { // Gly - no unit
         goto quitrout;
      }
      else if(tordata.sidedata[iside][0]==1){ // Ala and Pro - just one unit
         for(i=0;i<=tordata.sidedata[iside][3];i++){
           count=count+1;      
         }    
      } //Ala and Pro
      else{ // not Ala
         for(i=ns+1;i<=ns+tordata.sidedata[iside][0]-2;i++){
            count=count+1;
         }
         i=ns+tordata.sidedata[iside][0]-1;
         for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
             count=count+1;
         }
      } // not Ala
      quitrout:;
      sidecount[rottype] = count;
//      printf("For rottype = %d  %s, atomcount = %d\n",
//             rottype,nameres[rottype],count);
   } // rottype
}  
void setsideptspdb(struct geometry& geo,struct griddata& den1,
        struct pdbfile& pdb,int ***set,int segn,int resn)
/*
Similar to setside pts -  but only side chain
*/
{
   int i,ii,is,val,maxbr,maxpts,ib;

   val = geo.setval;
   is = segn;
   i = resn;
         maxbr = pdb.resdat[is+1][i][0]; // this = nbrn+1
         if(maxbr>1){ // then there are side chains
#if 0 // Old just set first 3 pts infirst branch
            maxpts = pdb.resdat[is+1][i][1]; // # of pts on first side branch, including Ca
            if(maxpts>3)
               maxpts =3; // set to max of 3 (or 2 side atoms, not counting Ca)
            for(ii=1;ii<maxpts;ii++)
              setatom(den1,geo,set,pdb.sidech[is][i][1][ii],val);//side chain pos
#endif
#if 1  // set all side atoms
            for(ib=1;ib<maxbr;ib++){ // each branch
               maxpts = pdb.resdat[is+1][i][ib]; // # of pts on first side branch, including Ca
//printf("For res %d, branch = %d, maxpts = %d\n",i,ib,maxpts);
               for(ii=1;ii<maxpts;ii++) // do not need to set ii=0 pt, because starts from previous pt
                   setatom(den1,geo,set,pdb.sidech[is][i][1][ii],val);//side chain pos
            }
#endif
         } // if maxbr
}

void setptspdb(struct geometry& geo,struct griddata& den1,
        struct pdbfile& pdb,int ***set,int segn,int resn)
/*
Similar to setpts - but use pdb rather than geo -
Sets all grid points within geo.setrad of mainchain
   and some side chain atoms
Included side chains atom names vary for each residue
Include 3 atoms in first side branch (or max if max < 3)
  First atom in side is always CA- do not need again
  Also, O atom is sidech[ii][i][j=nbrn+1= maxbr][1]
In lev  1=ca,2=c,3=o;4=n,5=cb; 
In mainch[0,1,2] = N,Ca,C
*/
#if 0
Excerpt from drawsidechain:
         for(ii=1;ii<=pdb1.resdat[is+1][i][0];ii++) {
            if(ii==pdb1.resdat[is+1][i][0])/* C-O bond*/
                 glColor3f(redv[0],redv[1],redv[2]);
            else 
              glColor3f(pdb1.sidecolor[is][0],pdb1.sidecolor[is][1],pdb1.sidecolor[is][2]);
            glBegin(GL_LINE_STRIP);
            for(j=0;j<pdb1.resdat[is+1][i][ii];j++) {
                glVertex3f(pdb1.sidech[is][i][ii][j][0],pdb1.sidech[is][i][ii][j][1],pdb1.sidech[is][i][ii][j][2]);   
            }
#endif
{
   int i,ii,is,val,maxbr,maxpts,ib;

   val = geo.setval;
   is = segn;
   i = resn;
         for(ii=0;ii<=2;ii++) // Main chain atoms 
            setatom(den1,geo,set,pdb.mainch[is][i][ii],val);
         maxbr = pdb.resdat[is+1][i][0]; // this = nbrn+1
         setatom(den1,geo,set,pdb.sidech[is][i][maxbr][1],val);// this is O atoms
         if(maxbr>1){ // then there are side chains
#if 0 // Only set first 3 points on first branch
            maxpts = pdb.resdat[is+1][i][1]; // # of pts on first side branch, including Ca
            if(maxpts>3)
               maxpts =3; // set to max of 3 (or 2 side atoms, not counting Ca)
            for(ii=1;ii<maxpts;ii++)
              setatom(den1,geo,set,pdb.sidech[is][i][1][ii],val);//side chain pos
#endif 
#if 1  // set all side atoms
            for(ib=1;ib<maxbr;ib++){ // each branch
               maxpts = pdb.resdat[is+1][i][ib]; // # of pts on first side branch, including Ca
//printf("For res %d, branch = %d, maxpts = %d\n",i,ib,maxpts);
               for(ii=1;ii<maxpts;ii++) // do not need to set ii=0 pt, because starts from previous pt
                   setatom(den1,geo,set,pdb.sidech[is][i][1][ii],val);//side chain pos
            }
#endif
         } // if maxbr
}

void refineonepdbside(int writetofile,struct griddata& den1,struct geometry& geo,
     struct rotomer rot[],struct pdbfile& pdb,int segn,int resn)
/*
   Similar to buildseq: 
  Refines the pdbside chain for pdb[segn][resn].
  BUT: May calcualte about bad contacts with other residues so do may use set
    and set useset = 0 so that use that do not use set in refineonerot
NOTE:IF MOVE Cb - THIS PROCESS IS NON-LINEar AND DOES NOT CONVERGE TO  BEST ANSWER -
   But gives better answer each time called.
  because, initially, translate all side chains by delta Cb.
Next time called, starting with whole side chain rotated so
   corresponds to new Ca-Cb bond, so get different result.
If Call 4 or 5 times in a row - roughly converges to best value.
Note: First time called, try all rotomers, (set bestrotn = -1),
    subsequent call just use previous bestrotn  
*/
{
   int npts,nbrn,count,atomname,bestnbad,bestbadset,rottype,bestrotn;
   int count2;
//   int lev0,lev1
   int ifit,fitres,useset;
   int irepeat,nrepeat;
   float bestrot[NAA][3];  
   float bestqual,origminrigid;
   float minqual,impfr;
   int minbad,usenew,initnbad;
   int i,j,k,iseg,ires;
   float cabond,catocb[3],endpos[3],endpos2[3];
   int origoneside;

  useset = 1; // means do use set in refinement
//  useset = 0; // means do not use set in refinement 
 origoneside = geo.oneside;
 if(useset==1){ // set all pdb pts = -1 first time called
    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;
    }
    if(pdb.iset == 0){ // not yet set
       printf("SETTING POINTS FOR ALL RESIDUES IN PDB #1\n");
       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;
       geo.setval = -1;
       for(iseg=0;iseg<pdb.resdat[0][0][0];iseg++) // all segments
         for(ires=1;ires<=pdb.resdat[iseg+1][0][0];ires++) 
              setptspdb(geo,den1,pdb,geo.set,iseg,ires);
       pdb.iset = 1;
    }
 } // useset
 
  impfr = 1.1; // qual must improve by this much before use new sidechain
//impfr = 0.0; //For degugging - use all changes qual must improve by this much before use new sidechain
// nrepeat = 1;
// pdb.resdat[segn+1][resn][4] = 0; // New 19jun00  -- Now set in allocatelev()   0 - default - not changed FIXED - increase NBRN to 5 and change 100 to 4 WARNING - OVERFLOW - NOT USED ANYMORE indicates has been modidied
// minqual =  0.7*maidfile.denscale; // min allowed to save fit
 minqual =  0.2*maidfile.denscale; // min allowed to save fit
// minbad = 1; // mininmum allowed to save new fit
 minbad = 2; // mininmum allowed to save new fit
 usenew = 0; // indicates have not found good enough fit to use new
 if(useset==1){
      geo.setval = 1; // set side = 1 so that will not hit itself
      setsideptspdb(geo,den1,pdb,geo.set,segn,resn);  
 }
#if 0
   if(dopowell==1){
       tordata.maxaunits= 3;// small value becuause just use one res.
       tordata.maxsunits=5*tordata.maxaunits;
       allocatetordata(tordata);
   }
#endif
   geo.phenylden= 0; //NOTE: usually does not do any good to use = 1   if = 1 then (in avedenoneside) find density only over the last rigid ring of a.a. reisdue - e.ge. phenylaline
   origminrigid = geo.minrigid;
   geo.minrigid = 0.5*maidfile.denscale; //Try increasing from old value = 0;  the value used to define nbad
   geo.oneside = 1;  // Since use copyregfittogeo - run dynamics just on geo #1 side chain
   ifit = 1;
   bestrotn = -1;// first time called, try all rotomers
#if 1
   fitres = 2; // the fit # that resn is copied to.
   copyreglevtofit(fit[ifit],pdb,segn,resn-2,resn+2);  // copies resn to fit #2  
#endif      
#if 0
   fitres = 0; // the fit # that resn is copied to.
   copyreglevtofit(fit[ifit],pdb,segn,resn,resn);  // copies resn to fit #2        
#endif      
   rottype = -fit[ifit].nside[fitres];
   if( (rottype ==GLY)||(rottype ==ALA) ) 
      goto skiprefine;
if(useset!=1){
   if( (rottype==ARG)||(rottype==LYS) ){// save terminal atom pos
      nbrn = data[rottype-1][1];
      count=1;
      count2 = 0; 
      for(j=1;j<=nbrn;j++){ // find name of last atom
          count=count+1;      
          npts=data[rottype-1][count];
          for(k=1;k<=npts;k++) {
            count=count+1;
            atomname=data[rottype-1][count];
            if(atomname>5)
               count2 = count2+1;
          }
      }// j
      atomname = count2+5;
      for(k=0;k<=2;k++){
         endpos[k] = fit[ifit].pos[fitres][atomname][k];
//printf("fitres = %d atomname = %denpos[%d] = %5.3f\n",fitres,atomname,k,endpos[k]);
      }
   }
   geo.oneside = origoneside; // set back to default
}// useset

#if 1 // Find initial qual
float bestden,initqual;
int iside;
int movetofar;
      iside = geo.oneside-1;
      geo.dir = 0;// N to C direction
    copyregfittogeo(geo,fit[ifit],fitres-1,fitres+1);
//    copyregfittogeo(geo,fit[ifit],fitres,fitres);
      geo.begend[0]= geo.nend[0]+iside;
      geo.begend[1]=geo.begend[0]+2;
//      geo.begend[1]=geo.begend[0];
//printf("begend[0] = %d  begend[1] = %d\n",geo.begend[0],geo.begend[1]);
    ncopygeotoauto(geo,tordata); 
    avedenoneside(den1,geo,tordata,1,&bestden,&bestnbad); // find true density for side chain = 1 - 
    initqual = bestden - maidfile.denscale*geo.valueweight[0]*bestnbad; // Assume no bad contacts.
    initnbad = bestnbad;
//    printf("INITIAL: qual = %5.3f  bestden = %5.3f  bestnbad = %d\n",initqual,bestden,bestnbad); // bestden same as above
#endif
  if(geo.movecb==1)
     nrepeat = 3;
  else 
     nrepeat = 1;
   for(irepeat = 1;irepeat<= nrepeat;irepeat++){ // Call several times to converve - see above
     copyregfittogeo(geo,fit[ifit],fitres-1,fitres+1);
//     copyregfittogeo(geo,fit[ifit],fitres,fitres);
#if 1 // try all rotomers each time
     bestrotn = -1;// first time called, try all rotomers
#endif
     refineonerot(useset,den1,geo,rot,ifit,fitres,rottype,&bestqual,
         &bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
     copyreggeotofit(geo,fit[ifit],fitres-1);
     if(useset!=1){
       if( (rottype==ARG)||(rottype==LYS) ){// save terminal atom pos
        for(k=0;k<=2;k++)
           endpos2[k] = fit[ifit].pos[fitres][atomname][k];
        if(distv(endpos2,endpos)> 4.0){
           printf("End atom move  too far in refine = %5.3f\n",distv(endpos2,endpos));
           goto repeat;
        }
      }
     } // useset
   if(geo.movecb==1){
     for(k=0;k<=2;k++)
        catocb[k] = fit[ifit].pos[fitres][5][k]- fit[ifit].pos[fitres][1][k];
     cabond = sqrt(  sqr(catocb[0])+sqr(catocb[1])+sqr(catocb[2]) );
//printf("bestnbad = %d  cabond = %5.3f\n",bestnbad,cabond);
     if(cabond>=2.0)
        movetofar = 1;
   }  
   if(geo.movecb==1){
     if( (bestqual> minqual)&&(bestnbad<= minbad)&&(cabond < 2.0)&&(cabond> 1.0)&&
                 (bestqual > impfr*initqual) ) {
         minqual = bestqual;
         copygeotorotsav(geo,geo.oneside,bestrot); //geo.oneside = resn
         usenew = 1;
     }
   }
   else{// do not move cb
      if( (bestqual> minqual)&&(bestnbad<= minbad)&&(bestqual > impfr*initqual) ) {
         minqual = bestqual;
         copygeotorotsav(geo,geo.oneside,bestrot); //geo.oneside = resn
         usenew = 1;
     }
   }
   repeat:;
  } // irepeat
  if(usenew == 1){
    pdb.resdat[segn+1][resn][4] = 1; //1 - not changed FIXED - increase NBRN to 5 and change 100 to 4 WARNING - OVERFLOW - NOT USED ANYMORE indicates has been modidied
    copyrotsavtogeo(geo,geo.oneside,bestrot); 
    copyreggeotofit(geo,fit[ifit],fitres-1);
    copyonefittolev(fit[1],pdb,fitres,segn,resn);
#if 0 // check if final qual same as bestqual
    copyregfittogeo(geo,fit[ifit],fitres-1,fitres+1);
      geo.begend[0]= geo.nend[0]+iside;
      geo.begend[1]=geo.begend[0]+2;
printf("begend[0] = %d  begend[1] = %d\n",geo.begend[0],geo.begend[1]);
    ncopygeotoauto(geo,tordata); 
    avedenoneside(den1,geo,tordata,1,&bestden,&bestnbad); // find true density for side chain = 1 - 
     initqual = bestden - maidfile.denscale*0.5*bestnbad;
    printf("Final :qual = %5.3f  bestden = %5.3f  bestnbad = %d\n",initqual,bestden,bestnbad); // bestden same as above
#endif
  }
if(geo.movecb==1){
   for(k=0;k<=2;k++)
     catocb[k] = fit[ifit].pos[fitres][5][k]- fit[ifit].pos[fitres][1][k];
    cabond = sqrt(  sqr(catocb[0])+sqr(catocb[1])+sqr(catocb[2]) );
//  printf("cabond = %5.3f\n",cabond);
}
skiprefine:;
float frcha,densc;
FILE *fp;
    densc = maidfile.denscale;
    if(writetofile==1){ // save to file
       if(initqual>0)
          frcha = (bestqual-initqual)/initqual;
       else
          frcha = 100;
       fp =fopen("refine.data","a");
       if( (rottype ==GLY)||(rottype ==ALA) )
          fprintf(fp,"  %s %s  (NOT REFINED)\n",
                pdb.rname[segn][resn],pdb.resnumch[segn][resn]);
       else if(usenew!=1)
          fprintf(fp,"  %s %s  (NOT IMPROVED)\n",
                pdb.rname[segn][resn],pdb.resnumch[segn][resn]);          
       else{
          if(geo.movecb==1)
            fprintf(fp,"Res %s %s: %5.3f %5.3f %d %5.3f %d %5.3f\n",
                pdb.rname[segn][resn],pdb.resnumch[segn][resn],frcha,initqual/densc,initnbad,
                bestqual/densc,bestnbad,cabond);
          else{
//                              Percent change in qual  Init. qual  Init. nbad  Final qual Final nbad  Ca-Cb
               fprintf(fp," %s %s:     %5.3f               %5.3f        %d        %5.3f      %d\n",
                pdb.rname[segn][resn],pdb.resnumch[segn][resn],frcha,initqual/densc,initnbad,
                bestqual/densc,bestnbad);
          }
       }
       fclose(fp);
    } 
    geo.minrigid=origminrigid;
#if 0
    if(dopowell==1)
      deallocatetordata(tordata);
#endif
//    savefit();
  geo.oneside = 0;// default - run standard dynnaicms on all atoms (if = 5, then just dynamics on sidechain #5)
 if(useset==1){
      geo.setval = -1; // set side back to -1
      setsideptspdb(geo,den1,pdb,geo.set,segn,resn);  
 }
}

void refineallpdbside(struct griddata& den1,struct geometry& geo,
     struct rotomer rot[],struct pdbfile& pdb)
/* Refine all the side chains in pdb
  Write side chain refinment data to "refine.data" and the 
    the New pdb to "refine.pdb"
*/
{
    int i,segn,resn,origsetrad;
    float origweight[2];
    FILE *fp;

    fitnum[0]=1;
    fitnum[1] = 1; // renumber fits so run consecuatively
    allocatefit(fit[1],10);
    origsetrad = geo.setrad;
    geo.setrad = 2; // Use large value here because assume start with good fit
    fp = fopen("refine.data","w"); // erase old file and start new
    fprintf(fp,"%s\n","Data for each sidechain residue in PDB #1 that was refined\n");
    if(geo.movecb==1)
       fprintf(fp,"Fraction change in qual  Init. qual  Init. nbad  Final qual Final nbad  Ca-Cb bond change\n");
    else
       fprintf(fp,"Fraction change in qual  Init. qual  Init. nbad  Final qual Final nbad\n");
    fclose(fp);
    for(i=0;i<=1;i++)
       origweight[i] = geo.valueweight[i];
#if 1 //Use less weight for nbad for refining pdb
    geo.valueweight[0] = 0.1; // value = -(aveden - maidfile.denscale*(valueweight[0]*nbad + valueweight[1]*badset) );
    geo.valueweight[1] = 1.0;
#endif
    for(segn=0;segn<pdb.resdat[0][0][0];segn++){ // all segments
        fp = fopen("refine.data","a"); // erase old file and start new
        fprintf(fp,"SEGMENT # %d\n",segn+1);   
        fclose(fp);
        for(resn=1;resn<=pdb.resdat[segn+1][0][0];resn++) { 
//        for(resn=1;resn<=100;resn++) { 
            printf(" REFINING RESIDUE %d:\n",resn);
            refineonepdbside(1,den1, geo,rot,pdb,segn,
                   resn); // 1 means write
        } // resn
     }// segn
     printf(" For valueweight[0,1] = %5.3f %5.3f\n",geo.valueweight[0],geo.valueweight[1]);
     writelevtopdb(pdb,"refine.pdb");
     geo.setrad=origsetrad;
     printf(" DONE: Data written to file 'refine.data' \n");
     printf("       and New PDB written to file 'refine.pdb' \n");
     for(i=0;i<=1;i++)
       geo.valueweight[i]=origweight[i];
}

void initsetptsfit(struct griddata& den1,struct geometry& geo)
/* sets points for all existing fits
*/
{
   int i,j,k;
   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(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;
   geo.setval = -1;
   for(i=1;i<=fitnum[0];i++){ //Extend side and  set points for each fit
if(verbose>=2)
 fprintf(fpout,"Setting set[][][] for fit %d\n",i);
      sidenumregfit(1,fit[i].nres-2,den1,geo,fit[i]); // Extends out to max. of geo.maxside, if there is density
      if(geo.quit==1){
         if(verbose>=2)
           printf(" Resetting geo.quit = 0 after sidenumregfit() set it to 1\n");
         geo.quit=0;
      }
      setptsfit(den1,geo,fit[i],geo.set,2); //"2" set first 3 side  pts - "1" indicates set all side set = geo.setval = -1;
   }
}

void buildseq(struct griddata& den1,struct geometry& geo,struct rotomer rot[])
/*
NOTE: (I think this is OK now)  BE CAREFUL ABOUT USING 
   sidenumreg and setpts if fit too big to copy to geo - see testseq
NOTE:  For each residue, first set side pts value = 1, so does not hit itself
           then, resets value to -1 using the new AMINO ACID SIDE CHAIN ATOMS
           Thus, this will differe slightly from original set.
  Routine that finds best side chain, assuming the sequence
   assignment in known.
  Input: geo.seqnum =  total number of residues in sequence
         geo.seq[i] = residue type (e.g. ARG) of ith residue
         geo.resnum = sequence number of picked residue in fit
         geo.pickfitnum = number of fit
         geo.fnum = residue number of this fit
*/
{
   int ifit,resn,bestnbad,bestbadset,bestrotn,rottype;
   int n0,origmaxside;
   int useset = 1; // means use set when refining
   float bestqual,origminrigid;
   int origsetrad;

if(verbose>=2)
   printf("Running buildseq\n");
//   setuprotomer(rot); // Now called in initialize.c++
#if 0 // New - 19jun00 - now define globally
    valueweight[0] = 0.5; // value = -(aveden - maidfile.denscale*(valueweight[0]*nbad + valueweight[1]*badset) );
    valueweight[1] = 5.0;
#endif
   geo.phenylden= 0; //NOTE: usually does not do any good to use = 1   if = 1 then (in avedenoneside) find density only over the last rigid ring of a.a. reisdue - e.ge. phenylaline
   origmaxside = geo.maxside;
   geo.maxside=3;
   origsetrad = geo.setrad;
   geo.setrad=2; // In initialize geo.setrad= 1;
   initsetptsfit(den1,geo);
   origminrigid = geo.minrigid;
   geo.minrigid = 0.5*maidfile.denscale; //Try increasing from old value = 0;  the value used to define nbad
   ifit = geo.pickfitnum;
   n0= fit[ifit].nres-1; // residues run from 0,1,2, .. n0, side chains from 1,n0-1
    geo.oneside = 1;  // Since use copyregfittogeo - run dynamics just on geo #1 side chain
    for(resn=1;resn<=n0-1;resn++){// run through each residue with sidechains
         copyregfittogeo(geo,fit[ifit],resn-1,resn+1);
         rottype = geo.seq[geo.resnum + resn - geo.fnum];
//printf("For resn %d  - type = %d\n",resn,rottype);
         geo.setval = 1; // set side = 1 so that will not hit itself
         setoneside(geo.oneside+geo.nend[0],den1,geo,geo.set,2);
         bestrotn = -1; // This means try all rotomers.
         refineonerot(useset,den1,geo,rot,ifit,resn,rottype,&bestqual,
             &bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
         fit[ifit].irot[resn] = bestrotn;
//printf("For ifit = %d resn= %d  rottype = %d    bestnbad = %d bestqual = %5.3f\n",
//        ifit,resn,rottype,bestnbad,bestqual);
         copyreggeotofit(geo,fit[ifit],resn-1);
         geo.setval = -1;
         setoneside(geo.oneside+geo.nend[0],den1,geo,geo.set,2); //NOTE: this set will be slightly different from original set = -1 set back to -1
      } // resn
    geo.maxside=origmaxside;
    geo.setrad=origsetrad;
   geo.minrigid=origminrigid;
//  savefit(); // Save fit after each build so that use new set()
  geo.oneside = 0;// default - run standard dynnaicms on all atoms (if = 5, then just dynamics on sidechain #5)
  quitrout:;
}

#if 0 // Old version - uses arrays (qual,truetype, etc that may overflow
void testseq(struct griddata& den1,struct geometry& geo,struct rotomer rot[])
/*
Set debug =1 to write out resqual for all fits - for debugging
    call from fit() menu (drawrrotomer)
  Similar to buildseq, but for each fit res, routine tries all amino acids and 
    returns the bestqual for each residue and prints out.
*/
{
   int i,ifit,resn,bestnbad,bestbadset,bestrotn,rottype;
   int n0,origmaxside,origmovecb,nbad,totfitnum,origsetrad;
   float bestqual,origminrigid;
   float aveden,firstden,lastden;
   float qual[100][VAL+5][4];  // qual for each residue for each amino acid
   int truetype[100],ibad[100][VAL+5]; // the real amino acid for that residue
   int sidecount[VAL+5],ifit0,seq[20],totside;
   int useset = 1;
   FILE *fp,*fpform;
int iend,debug;
 
//   debug = 0;  //For case where want to adjust qual function set debug = 1, read all fits and assume sequence assigment is know
   debug = 1;  //For case where want to adjust qual function set debug = 1, read all fits and assume sequence assigment is know
origmovecb = geo.movecb;
geo.movecb=1; // NEW 15sep00 - allow Cb to move when testing all a.a
printf("Running testseq and writing output to 'resqual' geo.movecb = %d\n",geo.movecb);    
#if 0
   if(dopowell==1)
      allocatetordata(tordata);
#endif
   geo.phenylden= 0; //NOTE: usually does not do any good to use = 1   if = 1 then (in avedenoneside) find density only over the last rigid ring of a.a. reisdue - e.ge. phenylaline
   origmaxside = geo.maxside;
   geo.maxside=3;
   origsetrad = geo.setrad;
   initsetptsfit(den1,geo);
   origminrigid = geo.minrigid;
   geo.minrigid = 0.5*maidfile.denscale; //Try increasing from old value = 0;  the value used to define nbad
   countside(geo,tordata,rot,sidecount);  // counts the number of side chain atoms in each amino acid tpe

seq[0]=GLY;seq[1]=ALA;seq[2]=SER;seq[3]=CYS;seq[4]=THR;seq[5]=VAL;seq[6]=PRO;seq[7]=LEU;
seq[8]=ASP;seq[9]=ASN;seq[10]=ILE;seq[11]=MET;seq[12]=GLN;seq[13]=GLU;seq[14]=HIS;
seq[15]=PHE;seq[16]=TYR;seq[17]=LYS;seq[18]=ARG;seq[19]=TRP;

fp=fopen(maidfile.resqual,"w");

if(debug==1){
   fpform=fopen("resqualform","w");
   fprintf(fpform," For each res:  qual  nbad  aveden firstden lastden\n");
   ifit0 = 1;
   iend = fitnum[0];
   totfitnum = fitnum[0];
#if 1 // just run on one fit
ifit0 = iend = 81;
totfitnum = 1;
#endif
}
else{
  ifit0 = iend = geo.pickfitnum;
  totfitnum = 1;
}

fprintf(fp,"%d\n",totfitnum); // totfitnum = number of fits written to resqul
for(ifit=ifit0;ifit<=iend;ifit++){ // each fit
   printf(" Assigning value for fit %d\n",ifit);
   n0= fit[ifit].nres-1; // residues run from 0,1,2, .. n0, side chains from 1,n0-1
   if(n0>98){ 
     printf("ERROR: n0>98 in testseq\n");
     exit(0);
   }
//resn=17;
    geo.oneside = 1;  // Since use copyregfittogeo - run dynamics just on geo #1 side chain
    for(resn=1;resn<=n0-1;resn++){// run through each residue with sidechains
       if(verbose>=2)
          printf("Testing fits for resn = %d\n",resn);
       copyregfittogeo(geo,fit[ifit],resn-1,resn+1);
       geo.setval = 1; // set side = 1 so that will not hit itself
       setoneside(geo.oneside+geo.nend[0],den1,geo,geo.set,2);
       if(geo.val1 < 0){
          geo.resnum = -geo.val1;
          truetype[resn]= geo.seq[geo.resnum + resn - geo.fnum]; // actual amino acid type
          if( (truetype[resn]<ALA) || (truetype[resn]>VAL) )
            printf("WARNING - TRUETYPE =%d  HAS AN IMPOSSIBLE VALUE\n",truetype[resn]);
       }
       else 
          truetype[resn] = 0;
       for(rottype = ALA; rottype<= VAL; rottype++){
         copyregfittogeo(geo,fit[ifit],resn-1,resn+1); // copy each time because Cb may be moved
         bestrotn = -1; // this means try all rotomers
         refineonerot(useset,den1,geo,rot,ifit,resn,rottype,
               &bestqual,&bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
         if(1){ // find detail for all residues
            ncopygeotoauto(geo,tordata); 
            detaildenoneside(den1,geo,tordata,1,1,&aveden,&nbad,
                &firstden,&lastden,&totside); // find true density for side chain = 1 - 
            qual[resn][rottype][1] = aveden;
            qual[resn][rottype][2] = firstden;
            qual[resn][rottype][3] = lastden;
         }
         else{
            qual[resn][rottype][1] = -1;
            qual[resn][rottype][2] = -1;
            qual[resn][rottype][3] = -1;
         }
         qual[resn][rottype][0] = bestqual;
         ibad[resn][rottype] = bestnbad;
      }//rottype
      copyregfittogeo(geo,fit[ifit],resn-1,resn+1); // copy original fit to geo
      geo.setval = -1;
      setoneside(geo.oneside+geo.nend[0],den1,geo,geo.set,2); // set back to -1
    } // resn
    fprintf(fp,"%d %d\n",ifit,n0);
    for(resn=1;resn<=n0-1;resn++){// run through each residue with sidechains
       fprintf(fp,"%d %d\n",resn,truetype[resn]);
//       printf("resn = %d, true type = %d %s:\n",resn,truetype[resn],nameres[ truetype[resn] ]);
       for(rottype = ALA; rottype<=VAL; rottype++){
#if 1 // print out all details
          fprintf(fp,"%d %f %d %f %f %f\n",rottype,qual[resn][rottype][0],ibad[resn][rottype],
                qual[resn][rottype][1],qual[resn][rottype][2],qual[resn][rottype][3]);
#endif
#if 0 // just print qual and nbad
          fprintf(fp,"%d %f %d\n",rottype,qual[resn][rottype][0],ibad[resn][rottype]);
#endif
#if 0
          printf("   type = %d   %s:  qual = %5.3f nbad= %d aveden = %5.3f  first = %5.3f last = %5.3f\n",
            rottype,nameres[rottype],qual[resn][rottype][0],ibad[resn][rottype],
            qual[resn][rottype][1],qual[resn][rottype][2],qual[resn][rottype][3]);
#endif
       } // rottype
    }// resn
   if( (debug==1)&& (-fit[ifit].val1>0) ){ //Write information file for adjusting parameters  Only write formated output  if there is correct sequence value for it
    fprintf(fpform,"\n\n Fit number %d: %d sidechains\n",ifit,n0-1);
    for(resn=1;resn<=n0-1;resn++){// run through each residue with sidechains
       fprintf(fpform,"Residue #%d:  Truetype = %s\n",resn,nameres[ truetype[resn] ]);
       for(i=1;i<=19;i++){
         rottype = seq[i];
#if 1 // print out all details
          fprintf(fpform,"%s: %f %d %f %f %f\n",nameres[rottype],qual[resn][rottype][0],ibad[resn][rottype],
                qual[resn][rottype][1],qual[resn][rottype][2],qual[resn][rottype][3]);
#endif
#if 0 // just print qual and nbad
          fprintf(fpform,"%d %f %d\n",rottype,qual[resn][rottype][0],ibad[resn][rottype]);
#endif
       } // i
       if(ibad[resn][AA3]==0)
           fprintf(fpform,"GOOD RING DEN: last qual = %f \n",qual[resn][AA3][0]);
       else
           fprintf(fpform," BAD RING DEN\n");
    }// resn
   }//if debug= 1 and ( -fit[i].val1>0)
 } // ifit
 fclose(fp);
 if(debug==1) 
    fclose(fpform);   
 geo.maxside=origmaxside;
 geo.setrad=origsetrad;
 geo.minrigid=origminrigid;
 geo.oneside = 0;// default - run standard dynnaicms on all atoms (if = 5, then just dynamics on sidechain #5)
 geo.movecb = origmovecb;
}
#endif

#if 1 // NEW (10oct00) version - does not use arrays (qual,truetype, etc that may overflow
void testseq(struct griddata& den1,struct geometry& geo,struct rotomer rot[])
/*
If geo.movecb = 0, do not move cb
              = 1  move for all residues

Set debug =1 to write out resqual for all fits - for debugging
    call from fit() menu (drawrrotomer)
  Similar to buildseq, but for each fit res, routine tries all amino acids and 
    returns the bestqual for each residue and prints out.
*/
{
   int i,ifit,resn,bestnbad,bestbadset,bestrotn,rottype;
   int n0,origmaxside,origmovecb,nbad,totfitnum,origsetrad,origringang;
   float bestqual,origminrigid,origmaxcbmove;
   float aveden,firstden,lastden;
   float qual[VAL+5][4];  // qual for each residue for each amino acid
   int truetype,ibad[VAL+5]; // the real amino acid for that residue
   int sidecount[VAL+5],ifit0,seq[20],totside;
   int useset = 1;
   FILE *fp,*fpform;
int iend,debug;
 
   debug = 0;  //For case where want to adjust qual function set debug = 1, read all fits and assume sequence assigment is know
#if 0
   printf("WARNING WARNING WARNING!!!  debug set =1 for debugging in testseq\n");
   debug = 1;  //For case where want to adjust qual function set debug = 1, read all fits and assume sequence assigment is know
#endif
origmovecb = geo.movecb;
origmaxcbmove = geo.maxcbmove;
geo.maxcbmove= 0.5; // New 30oct00 - try increasing allowed movement
geo.movecb=1; // NEW 15sep00 - allow Cb to move when testing all a.a
//geo.movecb=0; // New  default do notmove - but move for ring compounds in findqual- allow movement only for rings New 8dec00 - do not move cb - seems to make shadb a lot worse
origringang = geo.ringang; //if = 0 only use the default ring angle, else try others
geo.ringang=1; // In this test routine - try 3 ring angles -
if(verbose>=1)
   printf("Running testseq (tordata.fixend = %d) and writing output to 'resqual' geo.movecb = %d and maxcbmove = %5.3f\n",
          tordata.fixend,geo.movecb,geo.maxcbmove);    
#if 0
   if(dopowell==1)
      allocatetordata(tordata);
#endif
   geo.phenylden= 0; //NOTE: usually does not do any good to use = 1   if = 1 then (in avedenoneside) find density only over the last rigid ring of a.a. reisdue - e.ge. phenylaline
   origmaxside = geo.maxside;
   geo.maxside=3;
   origsetrad = geo.setrad;
   initsetptsfit(den1,geo);
   origminrigid = geo.minrigid;
   geo.minrigid = 0.5*maidfile.denscale; //Try increasing from old value = 0;  the value used to define nbad
   countside(geo,tordata,rot,sidecount);  // counts the number of side chain atoms in each amino acid tpe

seq[0]=GLY;seq[1]=ALA;seq[2]=SER;seq[3]=CYS;seq[4]=THR;seq[5]=VAL;seq[6]=PRO;seq[7]=LEU;
seq[8]=ASP;seq[9]=ASN;seq[10]=ILE;seq[11]=MET;seq[12]=GLN;seq[13]=GLU;seq[14]=HIS;
seq[15]=PHE;seq[16]=TYR;seq[17]=LYS;seq[18]=ARG;seq[19]=TRP;

fp=fopen(maidfile.resqual,"w");

if(debug==1){
   fpform=fopen("resqualform","w");
   fprintf(fpform," For each res:  qual  nbad  aveden firstden lastden lastbad\n");
   ifit0 = 1;
   iend = fitnum[0];
   totfitnum = 0;
   for(ifit=ifit0;ifit<=iend;ifit++){// count fits with sequence assigned
     if(fit[ifit].val1<0)
       totfitnum = totfitnum+1;
   }
#if 0// just run on one fit
ifit0 = iend = 15;
totfitnum = 1;
#endif
}
else{
  ifit0 = iend = geo.pickfitnum;
  totfitnum = 1;
}

fprintf(fp,"%d\n",totfitnum); // totfitnum = number of fits written to resqul
for(ifit=ifit0;ifit<=iend;ifit++){ // each fit
   if( (debug==1)&&(fit[ifit].val1>=0) )
        goto nextifit;
   printf(" Assigning value for fit %d\n",ifit);
   n0= fit[ifit].nres-1; // residues run from 0,1,2, .. n0, side chains from 1,n0-1
   geo.oneside = 1;  // Since use copyregfittogeo - run dynamics just on geo #1 side chain
   if( (debug==1)&& (-fit[ifit].val1>0) ) //Write information file for adjusting parameters  Only write formated output  if there is correct sequence value for it
       fprintf(fpform,"\n\n Fit number %d: %d sidechains\n",ifit,n0-1);
   fprintf(fp,"%d %d\n",ifit,n0);
//for(resn=n0-1;resn<=n0-1;resn++){// run through each residue with sidechains
   for(resn=1;resn<=n0-1;resn++){// run through each residue with sidechains
       if(verbose>=2)
          printf("Testing fits for resn = %d\n",resn);
       copyregfittogeo(geo,fit[ifit],resn-1,resn+1);
       geo.setval = 1; // set side = 1 so that will not hit itself
       setoneside(geo.oneside+geo.nend[0],den1,geo,geo.set,2);
       if(geo.val1 < 0){
          geo.resnum = -geo.val1;
          truetype= geo.seq[geo.resnum + resn - geo.fnum]; // actual amino acid type
          if( (truetype<ALA) || (truetype>VAL) )
            printf("WARNING - TRUETYPE =%d  HAS AN IMPOSSIBLE VALUE\n",truetype);
       }
       else 
          truetype = 0;
       for(rottype = ALA; rottype<= VAL; rottype++){
//for(rottype = ILE; rottype<= ILE; rottype++){
#if 0 // Try allowing cb to move just for rings
         if( (rottype==PHE)||(rottype==TYR)||(rottype==HIS)||(rottype==TRP) )
            geo.movecb = 1; // allow cb to move so ring fits perfectly
         else
            geo.movecb=0; 
#endif
         copyregfittogeo(geo,fit[ifit],resn-1,resn+1); // copy each time because Cb may be moved
         bestrotn = -1; // this means try all rotomers
         refineonerot(useset,den1,geo,rot,ifit,resn,rottype,
               &bestqual,&bestnbad,&bestbadset,&bestrotn);// returns geo = best rotomer
         if(1){ // find detail for all residues
            ncopygeotoauto(geo,tordata); 
#if 0 //Used before 18jan01 - does not check for bad contacts
            detaildenoneside(den1,geo,tordata,1,1,&aveden,&nbad,
                &firstden,&lastden,&totside); // find true density for side chain = 1 - 
#endif
#if 1 //Used after 18jan01 -  checks for bad contacts
            detaildenonesideset(den1,geo,tordata,1,1,&aveden,&nbad,
                &firstden,&lastden,&totside,1); //check for bad contact for i>=1  find true density for side chain = 1 - 
#endif
            qual[rottype][1] = aveden;
            qual[rottype][2] = firstden;
            qual[rottype][3] = lastden;
         }
#if 0 // This is no longer used with if(1)
         else{
            qual[rottype][1] = -1;
            qual[rottype][2] = -1;
            qual[rottype][3] = -1;
         }
#endif
         qual[rottype][0] = bestqual;
         if(rottype==PRO){ // use output of detaildenoneside
             ibad[rottype] = nbad; //8dec00 - these should be identical except for Proline case
             qual[rottype][0]=aveden - maidfile.denscale*geo.valueweight[0]*nbad ;  //Do not consider badset for these(alway bad for PRO)
         }
         else if(rottype==ILE){ // use output of detaildenoneside for nbad
             qual[rottype][0] = bestqual;
             ibad[rottype] = nbad; 
         }
         else if( (rottype==LYS)||(rottype==ARG) ){ // for lys and arg - if nbad in last unit - subtract 1 from nbad
              if(geo.res0[1]>0){ // at least 1 nbad in last unit
                 qual[rottype][0]= bestqual +  maidfile.denscale*geo.valueweight[0]; // equivalent of 1 bad
                 ibad[rottype] = bestnbad -1;
              }
              else {
                 qual[rottype][0] = bestqual;
                 ibad[rottype] = bestnbad;
              }
         }
         else{
           qual[rottype][0] = bestqual;
           ibad[rottype] = bestnbad;
         }
      }//rottype
      copyregfittogeo(geo,fit[ifit],resn-1,resn+1); // copy original fit to geo
      geo.setval = -1;
      setoneside(geo.oneside+geo.nend[0],den1,geo,geo.set,2); // set back to -1
      fprintf(fp,"%d %d\n",resn,truetype);
      for(rottype = ALA; rottype<=VAL; rottype++){
#if 1 // print out all details
          fprintf(fp,"%d %f %d %f %f %f\n",rottype,qual[rottype][0],ibad[rottype],
                qual[rottype][1],qual[rottype][2],qual[rottype][3]);
#endif
#if 0 // just print qual and nbad
          fprintf(fp,"%d %f %d\n",rottype,qual[rottype][0],ibad[rottype]);
#endif
       } // rottype
       if( (debug==1)&& (-fit[ifit].val1>0) ){ //Write information file for adjusting parameters  Only write formated output  if there is correct sequence value for it
          fprintf(fpform,"Residue #%d:  Truetype = %s\n",resn,nameres[ truetype ]);
          for(i=1;i<=19;i++){
             rottype = seq[i];
#if 1 // print out all details
              fprintf(fpform,"%s: %f %d %f %f %f %d\n",nameres[rottype],qual[rottype][0],ibad[rottype],
                   qual[rottype][1],qual[rottype][2],qual[rottype][3],geo.res0[1]); // res0[1] is nbad in last unit
#endif
#if 0 // just print qual and nbad
              fprintf(fpform,"%d %f %d\n",rottype,qual[rottype][0],ibad[rottype]);
#endif
          } // i
       }// debug==1
    } // resn
    nextifit:;
 } // ifit
 fclose(fp);
 if(debug==1) 
    fclose(fpform);   
 geo.maxside=origmaxside;
 geo.setrad=origsetrad;
 geo.maxcbmove=origmaxcbmove;
 geo.minrigid=origminrigid;
 geo.oneside = 0;// default - run standard dynnaicms on all atoms (if = 5, then just dynamics on sidechain #5)
 geo.movecb = origmovecb;
 geo.ringang=origringang; //if = 0 only use the default ring angle, else try others
}
#endif

#if 1 // New version - Use adjusted Ca-Ca distance
void lowerconnectnumca(struct geometry& geo,int ibestden,int *cnum,
       int *maxngap)
/*
  NEW return value of allowed maxngap.
 Moves along edited mainbone (= bone #1) and counts the number 
  of residues until reach last bone pt (=ext.length) = nterm Ca atom
  return number = cnum
Procedure: starting at pt0, move out until reach pt1
  where distance pt1 to pos0 >= geo.cadist.
Then find pos (= pos1) on the pt1-1 to pt1 branch where
   pos1 -pos0 = geo.cadist.
Then set pos0 = pos1 and repeat
After find pt where cdist > cadist
   then interpolate to find new pos0,
   And start next search at pos2 which is alread at 
      cdist = xlen
Note - since bone has been extended on next fit,
   there are bones pts for values > ext.length.
*/
{
   int pt0,k;
   float pos0[3],pos1[3],cdist,rvec[3],oldpos[3],oldcdist,fract; 

   pt0 = 0; // pt0 is the first Ca position
   for(k=0;k<=2;k++) // pos0 is the first Ca position
     pos0[k] = pos1[k] = ext.mainone.mainpt[pt0].pos[k];
   *cnum = 0;
   while(pt0< ext.length){
      cdist = 0.0;
      while(cdist<geo.cadist){
         pt0 = pt0+1;
         if(pt0>ext.length){
//            if(cdist> geo.cadist/3.0) // Use 1/3.0 to favor longer fits - since this usual underestimates
            if(cdist> geo.cadist/2.0)
               *cnum = *cnum +1;
printf("last addition: cdist = %5.3f, geo.cadist/3.0 = %5.3f\n",cdist,geo.cadist/3.0);
             goto atend;
         }              
         for(k=0;k<=2;k++){
           oldpos[k] = pos1[k];
           pos1[k] = ext.mainone.mainpt[pt0].pos[k];
         }
         oldcdist = cdist;
         cdist = distv(pos1,pos0);
      } 
      for(k=0;k<=2;k++)// at this point - have reach pos with dist > cadist -interpolate tofind caa pos
         rvec[k] = pos1[k]-oldpos[k];
#if 0 // OLD before 4oct00 - should not normalize - should use absolute rvec
      mag = sqrt( dot(rvec,rvec) );
      for(k=0;k<=2;k++)
         rvec[k] = rvec[k]/mag;
#endif
      fract = (geo.cadist - oldcdist)/(cdist - oldcdist);
      for(k=0;k<=2;k++) // the new capos - interpolated
         pos0[k] = oldpos[k]+fract*rvec[k];
      *cnum = *cnum +1; // coount ca's
      atend:; // jump to here is reach last point in above while(cdist..)
   } // while
#if 1 //Resolutin dependence - just crude guess - Add correction because bones tend to underestimate real distance - make resolution dependent
   if(geo.deldist> 1.4)
       *cnum = *cnum + *cnum/5;
   else if(geo.deldist > 1.3)
       *cnum = *cnum + *cnum/8;
   else if(geo.deldist > 1.2)
       *cnum = *cnum + *cnum/10;
   else if(geo.deldist > 1.1)
       *cnum = *cnum + *cnum/12;
   else if(geo.deldist > 1.0)
       *cnum = *cnum + *cnum/15;
#endif
   if(ibestden <= 1) // 0.8 or 1.0
       *maxngap = geo.maxngap;
   else if(ibestden ==2) // 1.2
       *maxngap = 1.5*geo.maxngap;
   else if (ibestden >2) // 1.4 or 1.6
       *maxngap = 2.0 *geo.maxngap;   
 }
#endif

void connectnumass(struct geometry& geo,struct griddata& den1,
         int fitext,int fitmet,int connum[],int ranki[],
         int ranki2[],int ranki3[],float rankf[],int ibestden,
         int cnum1,int nterm)
/*
 THIS IS CALLED ONLY IF ONE OF THE FITS HAS BEEN ASSIGNED!
fitext = nterm fit (if nterm = 0 then extend fitext to fitmet)
  fitmet = ctermfit  (if nterm = 1 then extend fitmet back to fitext)
  New version - only applied to case where one fit already assigned
  nterm = 0 -- cterm extension
  nterm = 1    nterm extension
*/
{
   float tvalue;
   int i,j,i1,tranki,tranki2,ii,seq0;
   int tnum1,mnum1,maxnum1;
   int allocate[5],nfits[MAXNF],seq1;
   int maxprint;

   ranki[0] = 0; // Number of saved ranked sequences
   connum[0] = 0;  // indicates fit not good enough to assign a definite gap
   allocate[0]=0;
   mnum1 = cnum1-geo.gapnum;
   if(mnum1<0)
      mnum1 = 0;
   maxnum1 = cnum1+geo.gapnum;
   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);
   }
   if(fit[fitext].val1 >=0){   
      geo.pickfitnum =fitext; // fit number
//      geo.resnum = geo.seq[0];//NEW - geo.resnum is is assigned value in testseq - true  give it the first value (NOT KNOWN IN THIS CASE) sequence  number for fit 1 res 1
      geo.fnum = 1; // res. number in fit
      testseq(den1,geo,rot);
      readresqual(maidfile.resqual,nfits); // only 1 fit is allocatere
//      printf("allocateing fitqual %d fits: = %d\n",nfits[0],nfits[1]);
      allocate[0]=allocate[0]+1; // counts fit allocated
      allocate[ allocate[0] ] = nfits[1];
//exit(0);
   }
   if(fit[fitmet].val1 >=0){    
      geo.pickfitnum = fitmet; // fit number
//      geo.resnum = geo.seq[0];//Dummy value  true  give it the first value (NOT KNOWN IN THIS CASE) sequence  number for fit 1 res 1
      geo.fnum = 1; // res. number in fit
      testseq(den1,geo,rot);
      readresqual(maidfile.resqual,nfits); 
//      printf("allocateing fitqual %d fits: = %d\n",nfits[0],nfits[1]);
      allocate[0]=allocate[0]+1; // counts fit allocated
      allocate[ allocate[0] ] = nfits[1];
   }
   printf(" Connect (gap) num for  gap= %d \n",cnum1);
   if( (fit[fitext].val1<0)&&(fit[fitmet].val1<0) ){ //both fitext and fitmet has sequence assignment - allow larger gap
     printf("Both fits have sequence assigned - testing gap using gaprange = %d:\n",geo.gapnum+1);
     maxnum1 = maxnum1+1;  
     printf("    fitext (%d)= %d and fitmet (%d)=%d\n",
             fitext,-fit[fitext].val1,fitmet,-fit[fitmet].val1);
     connum[1] = -fit[fitmet].val1 +fit[fitext].val1 - fit[fitext].nres +1;
     if( (connum[1]>=mnum1) && (connum[1]<=maxnum1) ){ // connum consistent with cnum1
        printf(" Gap number  = %d for assigment is consistent with bones\n",connum[1]);
        ranki2[1]=connum[1]; // the ideal gap
        connum[0]=1;
     }
     else {
       printf(" Gap number  = %d for assigment is NOT consistent with bones\n",connum[1]);
       connum[0] = -300;
     }
     goto quitrout; // Since ranki not define - skip printf("   seq. # = %d  gapnum = %d below  
  }
  else if(nterm==0){ // cterm extension
    if( fit[fitmet].val1>=0 ){ // this should always be case
      printf(" Nterm fit has already been assigned sequence# = %d\n",-fit[fitext].val1);
      if(verbose>=1)
         printf(" Test gap for cterm extension: range %d to %d\n",mnum1,maxnum1);
      seq0 = -fit[fitext].val1;
      for(tnum1 = mnum1; tnum1 <= maxnum1;tnum1++){
         seq1 =seq0+fit[fitext].nres+tnum1-1;
         if( (seq1+fit[fitmet].nres-3)>geo.maxseq){ // seq greater than max res. number in sequcne - 
           if(verbose>=1)
              printf("BAD test gap because seq1+fit[fitmet].nres-3=%d > geo.maxseq=%d\n",seq1+fit[fitmet].nres-3,geo.maxseq);
           goto nexttnum;
         }
         findqual(geo,fitmet,seq1,&tvalue);  // add to input tvalue
         tvalue = -tvalue; // Want to find smallest tvalue
         rankqual(MAXRANK,ranki,ranki2,rankf,seq1,tvalue,tnum1);
         nexttnum:;
      } //tnum1
    } // val1>0
    printf("Nterm fit (%d)already assiqned seq = %d: RANKING FOR Cterm fit (%d)=  \n",
            fitext,-fit[fitext].val1,fitmet);
    if(ranki[0]>7)
        maxprint = 7;
    else 
        maxprint = ranki[0];
    for(ii=1;ii<=maxprint;ii++)
         printf("   seq. # = %d  gapnum = %d  tvalue = %5.3f\n",
             ranki[ii],ranki2[ii],-rankf[ii]);
     if( (-rankf[1]> 25.0)||(ranki[0]==0) ){
        if(verbose>=1)
           printf("Bad connection because top ranking is > 25\n");
        connum[0]= 0;
     }
#if 0
     else if(-rankf[1]> 10.0){ // 8mar01 - if best rank >10 find gaps with in 10 of best
      for(ii=1;ii<=ranki[0];ii++){ // find gaps within 10 of best
         if( -rankf[ii]< -rankf[1] + 9.9){
           connum[0]=connum[0]+1;
           connum[ connum[0] ] = ii;
         }
       }
     }
#endif
     else{
       for(ii=1;ii<=ranki[0];ii++){ // find gaps within 4 of best
         if( -rankf[ii]< -rankf[1] + 3.9){
           connum[0]=connum[0]+1;
           connum[ connum[0] ] = ii;
         }
       }//ii
       for(ii=1;ii<=ranki[0];ii++){ // find gaps with tvalue of <7 and add to above list
         if( -rankf[ii]< 7.0){
            for(j=1;j<=connum[0];j++)
              if(connum[j] == ii) // already in list
                  goto nextii;
         }
         else
            goto quitii;
         connum[0]=connum[0]+1; // if reach here, not on list, add to list
         nextii:;       
       }// ii
       quitii:;
     } // else
  }// nterm = 0
  else{ //nterm =1, nterm extension 
     printf(" Cterm fit has already been assigned sequence# = %d\n",-fit[fitmet].val1);
     for(tnum1 = mnum1; tnum1 <=maxnum1;tnum1++){
        seq0 = -fit[fitmet].val1 -fit[fitext].nres -tnum1 +1;
        if(seq0<geo.seq[0]){ // seq0 less than minimum res. number in sequcne - 
           if(verbose>=1)
              printf("BAD test gap because seq0=%d < geo.seq[0]=%d\n",seq0,geo.seq[0]);
           goto nexttnum2;
        }
        findqual(geo,fitext,seq0,&tvalue);  // add to input tvalue
        tvalue = -tvalue; // Want to find smallest tvalue
        rankqual(MAXRANK,ranki,ranki2,rankf,seq0,tvalue,tnum1);
        nexttnum2:;
     } // tnum
     printf("Cterm fit (%d)already assiqned seq = %d: RANKING FOR  fitmet (%d)  \n",
            fitmet,-fit[fitmet].val1,fitext);
    if(ranki[0]>7)
        maxprint = 7;
    else 
        maxprint = ranki[0];
     for(ii=1;ii<=maxprint;ii++)
         printf("   seq. # = %d  gapnum = %d  tvalue = %5.3f\n",
             ranki[ii],ranki2[ii],-rankf[ii]);
     if(-rankf[1]> 25.0){
        if(verbose>=1)
           printf("Bad connection because top ranking is > 25\n");
        connum[0]= 0;
     }
#if 0 
     else if(-rankf[1]> 10.0){ // 8mar01 - if best rank >10 find gaps with in 10 of best
      for(ii=1;ii<=ranki[0];ii++){ // find gaps within 10 of best
         if( -rankf[ii]< -rankf[1] + 9.9){
           connum[0]=connum[0]+1;
           connum[ connum[0] ] = ii;
         }
       }
     }
#endif
     else{
       for(ii=1;ii<=ranki[0];ii++){ // find gaps within 4 of best
         if( -rankf[ii]< -rankf[1] + 3.9){
           connum[0]=connum[0]+1;
           connum[ connum[0] ] = ii;
         }
       }
       for(ii=1;ii<=ranki[0];ii++){ // find gaps with tvalue of <7 and add to above list
         if( -rankf[ii]< 7.0){
            for(j=1;j<=connum[0];j++)
              if(connum[j] == ii) // already in list
                  goto nextii2;
         }
         else
            goto quitii2;
         connum[0]=connum[0]+1; // if reach here, not on list, add to list
         nextii2:;       
       }// ii
       quitii2:;
     }// else
  } // nterm = 1
  for(ii=1;ii<=connum[0];ii++){ //Rank in order of gap length Simple sort routine - make connum[0] passes, first pass put shortest in first place, second pass place next shortes in second place, etc.
     for(i1=ii+1;i1<=connum[0];i1++){
        if(ranki2[i1]<ranki2[ii]){ // then interchange
             tranki = ranki[ii];
             tranki2 = ranki2[ii];
             ranki2[ii]= ranki2[i1];
             ranki[ii] = ranki[i1];
             ranki2[i1]=tranki2;
             ranki[i1]=tranki;
        }// interchange
     }//i1
  }//ii
  if(connum[0]>0){
     if( (-rankf[1]>15.0) && (ranki2[1]> 10) ) {
       if(verbose>=1)
         printf("SETTING connum[0]=0 because rankf[1]=%5.3f<15 and gap =%d >10\n",rankf[1],ranki2[1]);
       connum[0]=0;
       goto quitrout;
     }
  }
  if(verbose>=1){
    printf("Output of connectnjum:\n");
    for(ii=1;ii<=connum[0];ii++)
         printf("   seq. # = %d  gapnum = %d \n",
             ranki[ii],ranki2[ii]);
  }
  quitrout:;
  for(i=1;i<=allocate[0];i++){
      freefitqual(fit[ allocate[i]]);
  }

}

void connectnum(struct geometry& geo,struct griddata& den1,
         int fitext,int fitmet,int fit3,int connum[],int ranki[],
         int ranki2[],int ranki3[],float rankf[],int ibestden,
         int cnum1,int cnum2,int nterm)
/*
NEW 28sep00 - THIS IS CALLED ONLY FOR CASE WHERE geo.fitn is not assigned -
 DELTETED OTHER CASES.(see assignseq.c28sep00)
NEW - 12oct00 - allow case where fitmet already assigned - an extension
    from one direction may work when other direction failed.
  ONLY TRY THE BEST FIT.
 If ext.usenfits==3 then try use 3 fits to try to assign sequence
Only called if fitext has val1>0 - i.e. not assigned.
Only consider case where fitext and fitmet are not assigned.
   (If fitmet is assigned then should have previously tried (and failed) nterm
      extension from fitmet to fitext.)
If fit3 is assigned, then require that gap3 be consistent with assignment.
connum[0] = (only) -1 bad, or 1 found unique fit or if fitmet assigned, may return more than one
NEW stronger Criteria for assigning unique gap:
  1) If more than 1 fit has quality < 7.0:
      a) If cnum<=4:(which implies that it is probably the correct gap)
           Require that best gap is ideal gap, and second best fit with ideal gap >geo.rankgap  of first 
      b)   If cnum>4, then this is ambibuous case, do not use  
  2) If only 0 or 1 has quality < 7.0:
      Use the best gap if it is <rankmax and more than rankgap greater than second best.

         ranki[1,2..connum[0] = starting sequence number
         ranki1[1,2..]  = gap number
  NEW: IF more than 1 ranking has tvalue < 7 - Then:
     If ideal gap <=3 and there is an ideal gap with tvalue < 7
        and no other ideal gap within 3 of this ideal gap -
        then return connum[0] = 1 and gap = ideal gap.
     If ideal gap >4 - then treat this as an ambigous case
        and return connum[0]=0
NEW: 11jul00- ibestden is the the bone density used to find edited main
       (4 = 1.6, 0 = 0.8)
Routine to find the number of residues between two fits connected by an edited main
  fitext is the fit extended from cterm end
  fitmet in the fit that is met (meet the n term end of fitmet)
     Residue number  1 to fitext.nres to fitext.nres+ngap = to fitext.nres +ngap + fitmet.nres
Procedure:
   NOTE - new - lowerconnectum call before entering routine
      First call lowerconnectnum() - this return cnum = miminmub possible number of
        connecting residues.  Assume the correct num is  cnum-1 < = cnum(1.5)
   Then check if sequence has been assigned for either fit
     (if fit.val1 <0 then sequence of fit.res[1] =  - fit.val1)
   If not assigned - then run testseg for it.
Note program is written as if more than one fit can be indivually ranked
   However it only uses 1.
*/
{
   float tvalue,tvalue1,tvalue2,tvalue3,tvalue0;
   int i,ii,j,i1,seq0,seq1,seq2,nfits[MAXNF];
   float rankmax;
   int tnum1,totnum1,mnum1,mnum2,maxnum1,maxnum2;
   int allocate[5],tranki,tranki2;
   int tnum2,maxprint,test3gap;
   int icount;
   float f2;

   test3gap = 0; //default - do not test if third gap is OK (used if third fit already assigned)
   ranki[0] = 0; // Number of saved ranked sequences
   connum[0] = -1;  // indicates fit not good enough to assign a definite gap
   allocate[0]=0;
   if(ext.usenfits==3) // value is sum of 3 individual values (NOTE: rank values are normalized to length of 10?
     rankmax = 1.5*geo.rankmax;
   else // value is sum of 2 individual values
     rankmax = geo.rankmax;
   mnum1 = cnum1-geo.gapnum;
   if(mnum1<0)
      mnum1 = 0;
   maxnum1 = cnum1+geo.gapnum;
   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);
   }
   if( (ext.usenfits==3)&&(fit[fit3].val1<0) ) { //usenfits = 3 and  third fit is already assigned - rank first 2 fits
        if(verbose>=1)
           printf("Fit3 (%d) has seq = %d already assigned - see if compatible\n",fit3,-fit[fit3].val1);
        test3gap = 1; // Indicates should test third gap below
        ranki3[0] = -100; // returned value - indicates test3gap =1
   }
   if(fit[fitext].val1 >=0){ // this should always be case -sequence not assigned for this fit     
      geo.pickfitnum =fitext; // fit number
//      geo.resnum = geo.seq[0];//NEW - geo.resnum is is assigned value in testseq - true  give it the first value (NOT KNOWN IN THIS CASE) sequence  number for fit 1 res 1
      geo.fnum = 1; // res. number in fit
      testseq(den1,geo,rot);
      readresqual(maidfile.resqual,nfits); // only 1 fit is allocatere
//      printf("allocateing fitqual %d fits: = %d\n",nfits[0],nfits[1]);
      allocate[0]=allocate[0]+1; // counts fit allocated
      allocate[ allocate[0] ] = nfits[1];
//exit(0);
   }
   if(fit[fitmet].val1 >=0){ // this is not always the case - sequence not assigned for this fit     
      geo.pickfitnum = fitmet; // fit number
//      geo.resnum = geo.seq[0];//Dummy value  true  give it the first value (NOT KNOWN IN THIS CASE) sequence  number for fit 1 res 1
      geo.fnum = 1; // res. number in fit
      testseq(den1,geo,rot);
      readresqual(maidfile.resqual,nfits); 
//      printf("allocateing fitqual %d fits: = %d\n",nfits[0],nfits[1]);
      allocate[0]=allocate[0]+1; // counts fit allocated
      allocate[ allocate[0] ] = nfits[1];
   }
   if( (ext.usenfits==3)&&(fit[fitmet].val1 >=0)&&(fit[fitext].val1 >=0)&&(test3gap==0) ){// test fit3
      geo.pickfitnum = fit3; // fit number
      geo.fnum = 1; // res. number in fit
      testseq(den1,geo,rot);
      readresqual(maidfile.resqual,nfits); 
//      printf("allocateing fitqual %d fits: = %d\n",nfits[0],nfits[1]);
      allocate[0]=allocate[0]+1; // counts fit allocated
      allocate[ allocate[0] ] = nfits[1];
      mnum2 = cnum2-geo.gapnum;
      if(mnum2<0)
          mnum2 = 0;
      maxnum2 = cnum2+geo.gapnum;
   }
   printf(" Connect (gap) num for first gap= %d second gap =%d \n",cnum1,cnum2);
   if( (fit[fitext].val1>=0)&&(fit[fitmet].val1<0) ){ //NOTE THIS CASE SHOULD ONLY OCCUR FOR nterm = 0 (cterm extensio)See if gap compatible with assignment - This condition not really necessary because this connectin must have been already tried from other direction
      if(verbose>=1)
         printf("fitext not assigned and fitmet assigned seq = %d  nterm = %d\n",fit[fitmet].val1,nterm);
      for(tnum1 = mnum1; tnum1 <= maxnum1;tnum1++){
        if(nterm==0){ // cterm extension
          seq0 = -fit[fitmet].val1 -fit[fitext].nres -tnum1 +1;
          if(seq0<geo.seq[0]){ // seq0 less than minimum res. number in sequcne - 
             if(verbose>=1)
                printf("Not trying seq0=%d because less than allowed minmum = %d\n",seq0,geo.seq[0]);
             goto nexttnum;
          }
        }
        else{ // nterm extension  - THIS SHOULD NEVER HAPPEN!!
           seq0 = -fit[fitmet].val1+fit[fitmet].nres+tnum1-1;
           if( (seq0+fit[fitext].nres-3)>geo.maxseq){ // seq greater than max res. number in sequcne 
             if(verbose>=1)
                printf("Not trying seq0=%d because greater that allowed max = %d\n",seq0,geo.maxseq);
              goto nexttnum;
           }
        }
        findqual(geo,fitext,seq0,&tvalue);  // add to input tvalue
        tvalue = -tvalue; // Want to find smallest tvalue
        rankqual(MAXRANK,ranki,ranki2,rankf,seq0,tvalue,tnum1);
        nexttnum:;
     } // tnum
     printf("fitmet (%d)already assiqned seq = %d: RANKING FOR  fitext (%d)  \n",
            fitmet,-fit[fitmet].val1,fitext);
     if(ranki[0]>5)
        maxprint = 5;
     else 
        maxprint = ranki[0];
     for(ii=1;ii<=maxprint;ii++)
         printf("   seq. # = %d  gapnum = %d  tvalue = %5.3f\n",
             ranki[ii],ranki2[ii],-rankf[ii]);
     if( (-rankf[1]> 25.0)||(ranki[0]==0) ){
        if(verbose>=1)
           printf("Bad connection because top ranking is > 25\n");
        connum[0]= 0;
     }
     else{ //at least one good gap
       connum[0]=0;
       for(ii=1;ii<=ranki[0];ii++){ // find gaps within 4 of best
         if( -rankf[ii]< -rankf[1] + 3.9){
           connum[0]=connum[0]+1;
           connum[ connum[0] ] = ii;
         }
       }
       for(ii=1;ii<=ranki[0];ii++){ // find gaps with tvalue of <7 and add to above list
         if( -rankf[ii]< 7.0){
            for(j=1;j<=connum[0];j++)
              if(connum[j] == ii) // already in list
                  goto nextii;
         }
         else
            goto quitii;
         connum[0]=connum[0]+1; // if reach here, not on list, add to list
         nextii:;       
       }// ii
       quitii:;
     }// else
     for(ii=1;ii<=connum[0];ii++){ //Rank in order of gap length Simple sort routine - make connum[0] passes, first pass put shortest in first place, second pass place next shortes in second place, etc.
       for(i1=ii+1;i1<=connum[0];i1++){
        if(ranki2[i1]<ranki2[ii]){ // then interchange
             tranki = ranki[ii];
             tranki2 = ranki2[ii];
             ranki2[ii]= ranki2[i1];
             ranki[ii] = ranki[i1];
             ranki2[i1]=tranki2;
             ranki[i1]=tranki;
        }// interchange
       }//i1
     }//ii
   } // fitmet already assigned
   else if( (fit[fitext].val1>=0)&&(fit[fitmet].val1>=0) ){ //  Then try all sequence in pdbseqeucne
      if(verbose>=1)
         printf(" Neither of first 2 fits have sequence assigned\n");
//mnum1=maxnum1=1;  // for debugging
      for(tnum1 = mnum1; tnum1 <= maxnum1;tnum1++){
         totnum1 = fit[fitext].nres + fit[fitmet].nres - 1 +tnum1 +2; 
//geo.seq[0]=81; ;totnum1=geo.seq[0]-geo.seqnum;  // for debugging
         for(seq0=geo.seq[0];seq0<=geo.seqnum-totnum1;seq0++){ //shad starts with seq.# =12  pdb. seq. assigned to res 1
            findqual(geo,fitext,seq0,&tvalue1);
            seq1 =seq0+fit[fitext].nres+tnum1-1; 
            findqual(geo,fitmet,seq1,&tvalue2);  // add to input tvalue
            tvalue = -tvalue1-tvalue2; // Want to find smallest tvalue
            if( (ext.usenfits!=3)||(test3gap==1) ){
                  rankqual(MAXRANK,ranki,ranki2,rankf,seq0,tvalue,tnum1);
            }
            else{ // try 3 fits - the usual case, fit3 not assigned
                  tvalue0 = tvalue;
                  if( -tvalue0 < geo.rankmax){ //Only look at third fit if first 2 satisfy minimum criteria
                    for(tnum2 = mnum2; tnum2 <= maxnum2;tnum2++){
                       seq2 =seq1+fit[fitmet].nres+tnum2-1; 
                       if(seq2+fit[fit3].nres-3> geo.seqnum)
                          goto nextseq0;
                       else{ // possible assignment
                         findqual(geo,fit3,seq2,&tvalue3);  // add to input tvalue
                         tvalue = -tvalue1-tvalue2-tvalue3;
                         rankqual3(MAXRANK,ranki,ranki2,ranki3,rankf,seq0,tvalue,tnum1,tnum2);
                       }
                    } // tnum2
                 } // < rankmax
            }//try 3rd fit
            nextseq0:; // skip to here is sequence overlaps previous assigned fit
         } // seq0
      } //tnum1
      if(ranki[0]==0){
         printf("Could not assign sequence becaue ranki[0]=0!\n");
         goto quitrout;
      }
      else if(ranki[0]> 5)
         maxprint = 5;
      else 
         maxprint = ranki[0]; 
      if( (ext.usenfits !=3)||(test3gap==1) ){
         printf(" RANKING FOR fitext= %d fitmet =%d with totres = %d:\n",
            fitext,fitmet,fit[fitext].nres+fit[fitmet].nres-4);
         for(ii=1;ii<=maxprint;ii++)
            if(verbose>=1)
              printf("   seq. # = %d (%d)  gapnum = %d  tvalue = %5.3f\n",
                   ranki[ii],ranki[ii]+fit[fitext].nres+ranki2[ii]-1,ranki2[ii],-rankf[ii]);
      }
      else{ // print results for 3 fits
         printf(" RANKING FOR fitext= %d fitmet =%d fit3=%d with totres = %d:\n",
            fitext,fitmet,fit3,fit[fitext].nres+fit[fitmet].nres+fit[fit3].nres-6);
         for(ii=1;ii<=maxprint;ii++)
           if(verbose>=1)
            printf("   seq. # = %d (%d) (%d)  gap1 = %d gap2= %d tvalue = %5.3f\n",
             ranki[ii],ranki[ii]+fit[fitext].nres+ranki2[ii]-1,
             ranki[ii]+fit[fitext].nres+ranki2[ii]+fit[fitmet].nres+ranki3[ii]-2,
             ranki2[ii],ranki3[ii],-rankf[ii]);
      }
      if(test3gap==1){// Case where need 3 fits and third fit already assigned - best fit for first 2 must be compatible with third gap
          seq1 = ranki[1]+fit[fitext].nres+ranki2[1]-1; //best assignent for fit2 
          seq2 = -fit[fit3].val1;  // assigned seq. value for fit3
          tnum2 = seq2 - seq1 -fit[fitmet].nres +1;  // required third gap for best fit
          if( (tnum2>=mnum2)&& (tnum2<=maxnum2) ){ // if gap compatible with fit3 assignment
             if(verbose>=1)
               printf("Third fit already assigned and is Compatible with best fit\n"); 
             connum[1] = ranki2[1];
             connum[0]= 1;
          }
          else{
             if(verbose>=1)
               printf("Third fit is not Compatible with best fit\n"); 
          }              
      }
      else{ // test3gap = 0 - normal case
        icount =  0;
        for(ii=1;ii<=ranki[0];ii++){ // count value < 7 
         if( -rankf[ii]<= 7.0){
           icount = icount +1;
           if( -rankf[ii]>7)
              goto quitcount;
         }
        } // ii
        quitcount:;
        if(icount>=2){ // 2 or more values with very good rank
         if( (cnum1<=4)&&(ranki2[1]==cnum1) ){ // if short and best gap = ideal gap
            f2 = 1000.0;  // default value of second best with ideal gap
            for(ii=2;ii<=ranki[0];ii++){  // find next best with ideal gap
               if(ranki2[ii] == cnum1){ 
                   f2 = -rankf[ii];
                   goto quiti0count;
               }
            }
            quiti0count:;
            if(f2> -rankf[1]+geo.rankgap){ // then use best value
               connum[1] = ranki2[1];
               connum[0]= 1;
            }
         }
       } //icount
       else if( (-rankf[1]<=rankmax) && 
              (rankf[1]-rankf[2]) > geo.rankgap){ // If only 0 or 1 value < 7.0  then use best fit if <rankmax, and second best is much worse 
          connum[1] = ranki2[1]; // this is not needed??
          connum[0] = 1; // One good connectin
//        savefit();  // New - only save after assigning residues
       }
    }// test3gap = 0
  } // if neither fit assigned
  if(verbose>=1){
    if( (connum[0]!=1)&&(fit[fitmet].val1>=0) )// for this case, both fits unassigned, only allow conum[0]=1
       printf("Cannot assign fit because ranking ambiguous\n");
    else{
      printf("Output of connectnjum:\n");
      for(ii=1;ii<=connum[0];ii++)
         printf("   seq. # = %d  gapnum = %d \n",
             ranki[ii],ranki2[ii]);
    }
  }
  quitrout:;
  for(i=1;i<=allocate[0];i++){
      freefitqual(fit[ allocate[i]]);
  }
}


void assignone(struct geometry& geo,struct griddata& den1,int iexclude,int *nassign)
/*
  Just tries test one fit and sees if can assign sequence
  if iexclude = 1, then only test sequences that have not been
    assigned previously.
  This is only called after all posible sequences have been assigned and
     by using sym operators - assigned all possible unassigned
*/
{
   float tvalue,tvalue1;
   int ifit,i,ii,seq0,nfits[MAXNF];
   int *fitaa,count,ranki[MAXRANK+1],ranki2[MAXRANK+1]; //NOTE:BE SURE TO DELETE fitaa at end of routine
   float rankf[MAXRANK+1];
   int tnum1,totnum;
   int allocate[5];
   int maxprint;

   *nassign =0;
   printf("Trying to assign sequences to all unassgined fits using iexclude = %d,geo.rankmax = %5.3f,geo.rankgap= %5.3f\n",
                 iexclude,geo.rankmax,geo.rankgap);
   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);
   }
   tnum1 = 0;  // NOT used - just need dummy value for rankqual()
   if(iexclude==1){ //Find fits that have allready been assigned exclude seq0 based on other fit
     count = 0;
     for(i=1;i<=fitnum[0];i++) //first count then store fits that have sequence assigned
        if(fit[i].val1<0)
           count = count +1;
     fitaa = new int[count+4]; // allocate fitaa   
     fitaa[0] = 0;
     for(i=1;i<=fitnum[0];i++){ // store fits that have sequence assigned
        if(fit[i].val1<0){
           fitaa[0]= fitaa[0]+1;
           fitaa[ fitaa[0] ] = i;
        }
     } // i
   }// iexclude
   for(ifit=1;ifit<=fitnum[0];ifit++){
//for(ifit=11;ifit<=fitnum[0];ifit++){
      ranki[0] = 0; // Number of saved ranked sequences
      allocate[0]=0;
      if(fit[ifit].val1 >=0){ // sequence not assigned for this fit
        if(verbose>=1)     
           printf("Testing fit %d for assignment  using iexclude = %d\n",ifit,iexclude);
        geo.pickfitnum =ifit; // fit number
        geo.fnum = 1; // res. number in fit
        testseq(den1,geo,rot);
        readresqual(maidfile.resqual,nfits); // only 1 fit is allocatere
//      printf("allocateing fitqual %d fits: = %d\n",nfits[0],nfits[1]);
        allocate[0]=allocate[0]+1; // counts fit allocated
        allocate[ allocate[0] ] = nfits[1];
        totnum = fit[ifit].nres +2;
        rankf[1] = -10001.0;  // default value - indicates never ranked 
        for(seq0=geo.seq[0];seq0<=geo.seqnum-totnum;seq0++){ //shad starts with seq.# =12  pdb. seq. assigned to res 1
           if(iexclude==1){ //exclude seq0 based on other fit
             for(i=1;i<=fitaa[0];i++){ // do not test seq0 if it starts or end in a previous assigned sequence
               if( (seq0>= -fit[ fitaa[i] ].val1+1)&&
                      (seq0<= -fit[ fitaa[i] ].val1+ fit[fitaa[i] ].nres-3)) //starting residue overlaps assinged fit
                   goto nextseq0;
               if( (seq0+fit[ifit].nres-2 >= -fit[ fitaa[i] ].val1+1)&&
                      (seq0+fit[ifit].nres-2<= -fit[ fitaa[i] ].val1+ fit[fitaa[i] ].nres-3)) //last residue overlaps assinged fit
                   goto nextseq0;
             }
           }
           findqual(geo,ifit,seq0,&tvalue1);
           tvalue = -tvalue1; // Want to find smallest tvalue
           rankqual(MAXRANK,ranki,ranki2,rankf,seq0,tvalue,tnum1);
           nextseq0:;
        }// seq0
        if(verbose>=1){
          if(ranki[0]>5)
            maxprint = 5;
          else 
            maxprint = ranki[0];
          for(ii=1;ii<=maxprint;ii++)
           printf("   seq. # = %d    tvalue = %5.3f\n",
               ranki[ii],-rankf[ii]);
        }
        if( (-rankf[1]<=geo.rankmax) && 
              ( (rankf[1]-rankf[2]) >= geo.rankgap) &&
                 (-rankf[2]>10.0) ){ // Then assign sequence to this fit
           if(verbose>=1)
              printf("Assigning fit %d sequence number = %d\n",ifit,-ranki[1]);
           *nassign = *nassign+1;
           fit[ifit].val1 = -ranki[1];
           buildandrefineseq(den1,geo,rot,ifit);           
        }
        else if(rankf[1]< -10000){
          if(verbose>=1)
             printf("Sequence not assigned all seq0 already assigned\n");
        }
        else{
          if(verbose>=1)
            printf("Sequence cannot be assigned because best = %5.3f,nextbest = %5.3f\n",
                     rankf[1],rankf[2]);
        }
      } //val1>0
      for(i=1;i<=allocate[0];i++)
          freefitqual(fit[ allocate[i]]);
   }// ifit
   savefit();
}


void buildandrefineseq(struct griddata& den1,struct geometry& geo,
         struct rotomer rot[],int ifit)
/* For fit ifit - call buildseq and refine using torsional dynamics
   on all residues.
 NOTE: buildseq uses copyregfittogeo, so it can be used on very 
     long fits.  However, refinement assumes can refine all residues
     in geo (geo.fixend = 0) which means geo cannot be very long.
   Maybe should modify so that if too long, refine in segments.
*/
{
   float lastden;
   float origrfor,origtorconstend;
   int origfixend,origoneside;
   int nbad,quitcond,val[2];

    geo.pickfitnum = ifit; // fit number that is used in buildseq
    geo.resnum = -fit[geo.pickfitnum].val1 ; // pdb sequence number
    origfixend = geo.fixend;
    origoneside = geo.oneside;
    origrfor = tordata.rfor;
    origtorconstend = tordata.torconstend;
    val[0] =fit[geo.pickfitnum].val1; 
    val[1] =fit[geo.pickfitnum].val2; 
    if(geo.resnum <=0){// have not yet assigned seqence number
       printf("ERROR: Have not yest assigned seq. number in buildandrefine -skipping\n");
       goto quitrout;
    }
    else{ // build and refine
      geo.fnum = 1; // res. number in fit corresponding to above geo.resnum
      buildseq(den1,geo,rot); // This returns the new fit
      geo.fixend = 0; // run tordyn on all residues
      geo.oneside = 0;  // run on all sides
#if 1
if(verbose==2)
   savefit();
#endif
      copyfittogeo(geo,fit[geo.pickfitnum]); //NOTE ASSUMES FIT IS NOT TOO LONG      
      quitcond = 30; // run 30 rounds of refinement
      if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
               tordata.adata[1] = 2; 
 tordata.rfor = 0; // NEW 12dec00 - try using 0 random force here.
      tordynbat(quitcond,den1,geo,tordata,&lastden,&nbad);// "1" means use N term density as criteria for q;uiting
#if 0 // Use run round of dynamics with end density force increased
      tordata.torconstend = 10.0*origtorconstend;// MAY HELP IN SOME CASES - BUT NOT COMPLETELY TESTED New 21dec00 - try increasing density  force on end here
#endif
      tordynbat(quitcond,den1,geo,tordata,&lastden,&nbad);// "1" means use N term density as criteria for q;uiting
      printf(" Assigning residues to fit %d and refining using tordata.rfor/origrfor = %5.3f and torconstend/orig = %5.3f\n",
              ifit,tordata.rfor/origrfor,tordata.torconstend/origtorconstend);
      if(geo.quit==1){
         printf("NOTE - Not using torsion dyn because out of map range\n");
      }
      else
         copygeotofit(geo,fit[geo.pickfitnum],val); // Copy 
   }
   quitrout:;
    geo.fixend=origfixend;
    geo.oneside=origoneside;
    tordata.rfor=origrfor;
    tordata.torconstend=origtorconstend;
 }

#if 0 // Used for debugging - routine lowerconnectnum() has been deleted
void testgap(struct griddata& den1,struct geometry& geo)
/* routine to test assignment of ideal gap number
*/
{
   int i,ifit,idealgap,foundedit,ibestden,ibeg,iend,count;
   int isummary[50][7];
   int canum,cnum,maxngap;

   geo.fixend =2; // This means extend at c term end
   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);
    if(geo.fitn==0){
       ibeg= 1;
       iend = fitnum[0];
    }
    else{
       ibeg=iend = geo.fitn;
    }
    count=0;
    for(ifit = ibeg;ifit<=iend;ifit++){
       geo.fitn = ifit;
       batchextend(den1,geo,&foundedit,&ibestden);
       lowerconnectnumca(geo,ibestden,&canum,&maxngap); 
       if( (ext.length>0)&&(canum<=maxngap) ){
         if( (fit[geo.fitn].val1<0)&&(fit[ext.metend].val1<0) ){ //both fitext and fitmet has sequence assignment
            printf("Both fits have sequence assigned:\n");
            printf("    fitext (%d)= %d and fitmet (%d)=%d\n",
             geo.fitn,-fit[geo.fitn].val1,ext.metend,-fit[ext.metend].val1);
            idealgap = -fit[ext.metend].val1 +fit[geo.fitn].val1 - fit[geo.fitn].nres +1;
            printf(" gap number = %d  ideal gap = %d\n",cnum,idealgap);
         }
         count = count+1;
         isummary[count][1]=idealgap;
         isummary[count][2]=geo.fitn;
         isummary[count][3]=ext.metend;
         isummary[count][4]=cnum;
         isummary[count][5]=canum;
         isummary[count][6]=ibestden;
       }// if ext.length>0
    } 
    isummary[0][0]=count;
    for(i=1;i<=isummary[0][0];i++){
       printf("Fit %d met fit %d bestden=%d: my best gap = %d  ideal gap = %d oldgap=%d cagap = %d\n",
           isummary[i][2],isummary[i][3],isummary[i][6],isummary[i][1],isummary[i][4],isummary[i][5]);
    }
}
#endif

void findside_maindistfit(struct pepfit& fit,int ires,float *mindist)
/* Finds the minimum dist of residue side chain atoms from Cb2 on out
  Checks distance from Ca -C of previous and N- Ca-C of current residue and N-Ca of next residue
   fit.nside[fitloc] = -rottype;  
   geo.iside[geoloc] = -rottype; 
*/
{
   int i,ii,j,rottype,count,set[3][4],atomname,tres,useset,iatom;
   float dist;
  
   rottype = -fit.nside[ires];  
   if(rottype <=0) // not a residue
     goto quitrout;
   if( (rottype==ALA)||(rottype==VAL)||(rottype==PRO)||(rottype==ILE)||(rottype==LEU)
        ||(rottype==ASP)||(rottype==SER)||(rottype==THR)||(rottype==CYS)||(rottype==ASN) )
     goto quitrout;
   count=0; // count side atoms;  count = 1 corresponds to Cb2 
   for(ii=1;ii<=side[rottype][0][0]+1;ii++){ // copy each side chain branch
      for(j=0;j<side[rottype][ii][0];j++){
         atomname=side[rottype][ii][j+1];
         if(atomname>5){ // not part of mainchain or Cb
           count = count+1;
         }
      }
   }//ii
   set[0][0]=2;
   set[0][1]=1;//CA
   set[0][2]=2;//C
   set[1][0]=3;
   set[1][1]=1;//CA
   set[1][2]=2;//C
   set[1][3]=4;// N
   set[2][0]=2;
   set[2][1]=1;//CA
   set[2][2]=4;//N
   for(i=2;i<=count;i++){ //i=2 correspondns to Cb3
      for(tres = ires-1;tres<=ires+1;tres++){ // check dist for tres main chains
         useset = tres-ires+1;
         for(ii=1;ii<=set[useset][0];ii++){
            iatom = set[useset][ii];          
            dist = distv(fit.pos[tres][iatom],fit.pos[ires][count+5]);
            if(dist< *mindist){
               *mindist = dist;
printf(" mindist = %5.3f for res %d  rotype = %d count= %d\n",*mindist,ires,rottype,count);
            }//dist
         }//ii
      }//tres
   }//i
   quitrout:;
} 


void findside_maindisttor(struct geometry& geo,int ires,int iside,float *mindist)
/* 
Same as findside_maindistfit except use tordata for side and geo for main chain
ires = geo.resn number
iside usual = 1 and in is the sidenum in tordata
*/
{
   int i,ii,j,rottype,count,set[3][4],tres,useset,iatom;
   int ns;
   float dist;
  
   set[0][0]=2;
   set[0][1]=1;//CA
   set[0][2]=2;//C
   set[1][0]=3;
   set[1][1]=1;//CA
   set[1][2]=2;//C
   set[1][3]=4;// N
   set[2][0]=2;
   set[2][1]=1;//CA
   set[2][2]=4;//N
   rottype= tordata.sidedata[iside][2];
   ns=tordata.sidedata[iside][1];// start number for this chain = Ca
   if(rottype <=0) // not a residue
     goto quitrout;
   if( (rottype==ALA)||(rottype==VAL)||(rottype==PRO)||(rottype==ILE)||(rottype==LEU)
        ||(rottype==ASP)||(rottype==SER)||(rottype==THR)||(rottype==CYS)||(rottype==ASN) )
      goto quitrout;
   count=0; // count side atoms;  count = 1 corresponds to Cb 
   for(i=ns+1;i<=ns+tordata.sidedata[iside][0]-2;i++){ // atoms up to last unit
      count=count+1;
      if(count>=3){// Cb3 and greater
        for(tres = ires-1;tres<=ires+1;tres++){ // check dist for tres main chains
         useset = tres-ires+1;
         for(ii=1;ii<=set[useset][0];ii++){
            iatom = set[useset][ii];          
            dist = distv(geo.pos[tres][iatom],tordata.uside[i][0]);
            if(dist< *mindist){
               *mindist = dist;
//printf(" mindist = %5.3f for res %d  rotype = %d count= %d\n",*mindist,ires,rottype,count);
            }//dist
         }//ii
        }//tres
      }
   }
   i=ns+tordata.sidedata[iside][0]-1;
   for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
      count=count+1;
      if(count>=3){// Cb3 and greater
        for(tres = ires-1;tres<=ires+1;tres++){ // check dist for tres main chains
         useset = tres-ires+1;
         for(ii=1;ii<=set[useset][0];ii++){
            iatom = set[useset][ii];          
            dist = distv(geo.pos[tres][iatom],tordata.uside[i][j]);
            if(dist< *mindist){
               *mindist = dist;
//printf(" mindist = %5.3f for res %d  rotype = %d count= %d\n",*mindist,ires,rottype,count);
            }//dist
         }//ii
        }//tres
      }
   }
   quitrout:;
}

void testseqbatch(struct griddata& den1,struct geometry& geo)
/*Called from mainbatch -
  Runs  testseq() - which read fits in maidfit - and generates
     resqual file.
  then rankseq() - which uses resqual file to rank all the fits.
  NOTE: NEED TO HAVE SEQUENCE ASSIQNED IN VAL1 BEFORE RUNNING -
   (see rankside.notes for detail and  HOW TO ASSIGNE VAL1 AUTOMATICALLY
  Also prints resqualform which is easier to read. 
*/
{
     geo.pickfitnum = 2; // fit number - not used if debug = 1
     geo.resnum = 180;// sequence  number for fit 1 res 1  - not used if debug = 1
     geo.fnum = 1; // res. number in fit
     readfit(geo,maidfile.fitfile);
     if(geo.seqnum== -1){ // indicats pdb sequence not yet read - call readsequence
        if(readsequence(geo,maidfile.pdbsequence)!=1){ // indicates not successfully read
           printf("ERROR: one letter sequence file 'maidfile.pdbsequence' is missing\n");
           goto quitrout;
        }
     }
#if 1 // If debug =1 in testseq(), then writes out for all fits that have sequence assigned
     testseq(den1,geo,rot);
#endif
#if 1 // run rankseq after
     geo.seqnum = -1; // indicates pdbsequence not read
    rankseq(geo,maidfile.resqual);
#endif
quitrout:;
}
