/* tordyn.c++   Routines to run tordyn  - adapted from old "bonedyn" 
         directory that used the old maid routines

6mar00: Found a subtle error - did not set all initial
   velocities = 0 in tordynbat
20mar00:  Saved old version in tordyn.c20mar00
   Slightly modified ncopygeotoauto and 
       substituted ncopygeo.. for copygeo.. in tordynbat()
28JAN00: PROBLEM:
    IF tordata.sidedata[i][0] = 1  (ALA and PRO) then in runtordyn():
  for(i=1;i<=tordata.nunit1;i++){ 
      if(tordata.sidedata[i][0]>0){
         for(j=tordata.sidedata[i][1]+1;j<tordata.sidedata[i][1]+tordata.sidedata[i][0];j++){
  This code will do nothing if tordata.sidedata[i][0] = 1.
  For standard torsional dynamics, it would mean that Ala and Pro
     side chains would not be affect by local density.
 FIXED - SIMPLY SET  tordata.sidedata[i][0] = 2 for PRO and ALA
   have debugged for Ala, have not checked Pro, but it should work.
12nov99:  COMMENT CALL OF checkanglimit  IN tordynbat
8oct99:  Found major error in ncopygeotoauto for Cto N direction - cause problem in tordyn about pick:
   tordata.unit[0][i][0][k] should be tordata.unit[0][i][1][k]
   tordata.unit[0][i][0][k]=geo.pos[ii][5][k];  Cb pos 
  should be
   tordata.unit[0][i][1][k]=geo.pos[ii][5][k];  Cb pos 
  I do not know why this did not show up as problem when run tordyn at N term end.
  Also - I have not check if side exist before defining.
    and why problem just shows up now.
  copygeotoauto was OK so tordynbat runs OK.
9aug99   NEW  version -treat all terminal atoms as one rigid unit in ILE
    in copyonesidetogeo- neglects one degree of freedom
    ILE does not fit general rule that only unit with more than
    one atom is last unit.
    ---WORKS OK  IN TEST ON endodm.
On 17may99:  Added procedure for either running dynamics in middle of
    chain (geo.fixend = 3)  - and - same procedure can be used for "fusing" two chains.
Procedure: Dynamics around residue i:  (defined by sidechain = i)
    1) Set up to run dynanmics for Cterm end, from residue i-3 to i.
    2)  Run Nterm dyncamics for residue i+3 to i.
    3) NEW: Add force to both end residues atoms to minimize distance
           between atoms:  
         Cterm:  geo.nend[1]=i+1;
            Atoms Ca[i+1], N[i+1], C[i], O[i]
         Nterm:  geo.nend[0]=i-1;
            Atoms Ca[i-1], N[i-1], C[i], O[i]

Timing:
I) 1000 runs of copygeotoauto for Nterm dynamics
Start copygeotoauto  CPU time = 0
End copygeotoauto  CPU time = 90000
II) 1000 runs of runtordyn for Nterm dynamics
Start runtordyn  CPU time = 0
End tordyn  CPU time = 2570000
III) 1000 runs of runtordyn for Nterm dynamics
Start copyautotogeo  CPU time = 0
End copyautotogeo  CPU time = 10000
NOTE: copygeotoauto takes longer than autotogeo because must find
     angsign.
COPYING TAKES ABOUT 10/257 OF TIME OF TORDYN - 
    thus - About 10% slow down if copy back and forth between
    each call to tordyn when run in middle.
DEFINE
  float geo.fpos[2][7][3]; // pos of end residue - used when fuse two fits or run dyn. in center
  geo.fpos[0][][]  pos of N term atoms and fpos[1][][]  pos of C term atoms
    geo.fpos[0][Define  1=ntermca,2=c,3=o;4=n,5=cb 6=ctermca ]  (cb not used now)
USE SPECIAL ROUTINE FOR drawing geo when geo.fixend =3, draws both
   ends separately  

    
On may7,99 -TRIED:  Use cbdenreg to determine if cb in density and
     in runtordynnew - 
   change
       #if 0// leave out cbden
              natoms = 1;
       #endif
           natoms=2;
    to include cbden
 RESULT:  This maid trace worse - CHANGED BACK DID NOT USE

mar31,99 - int tordynbat -use initial position as initial bestden
March 5, 99:   Added following line before call to aventermden 
   (before - this was only used for avecterm)
   I really do not know if this cbden thing is working correctly
      for(i=0;i<NAFIT;i++)
         geo.cbden[i]=1;

Note on 23NOV98 - possible error in procedure -
  When finding D array (dd) use the last Y vector, which is found by
     integrating dY.  However, probably should go directly to new atoms positions,
     find dr and dw from atom postions (made using makealanine) and then
     use this dr and dw to find Y and D.
 However, it is not easy to find dw from atom positions.(dr is simple).
 ACTUALLY, I AM SURPRISED THAT MY PROCEDURE WORKS OK GIVEN THIS PROBLEM.
   I WOULD EXPECT THAT AFTER A LONG TIME, dY would not be consistent with
   actual velocities.
  However, when tested this by comparing actual c.m. velocities, with
    those estimated from Y, there was no evidence that error got 
    worse with time.

NOTE: New on 20Nov98  - modified frictional term so that is more accurate -
      keep track of velocity of each atom directly (in makeunitalanine)
      and use this velocity to find resistance.
New on 11 Nov 98  - Allow amino acid side chains:
 NOTE: THIS HAS NOT BEEN IMPLEMENTED FOR GENERAL CASE, JUST
         FOR CASE OF MOVING ONE SIDE CHAIN - SHOULD BE EASY TO GENERALIZE
  1) tordata.sunit now has dimension   sunit[resnum][10][3] to allow for tryptophane atom
  2) geo.iside = -rottype
  3) only the last unit has multiple atoms
  5) For each residue - (BRUTE FORCE APPROACH!!)
      have separate copygeotoauto and copyautotogeo:
    e.g.
       phenyalanine:
          tordata.uside[ns][0][k] =geo.pos[ii+1][5][k]; (Cb)
          tordata.uside[ns+1][0][k] =geo.pos[ii+1][6][k]; (CG)
          tordata.uside[ns+1][1][k] =geo.pos[ii+1][7][k]; (CD1)
          tordata.uside[ns+1][2][k] =geo.pos[ii+1][8][k]; (CE1)
          tordata.uside[ns+1][3][k] =geo.pos[ii+1][9][k]; (CZ)
          tordata.uside[ns+1][4][k] =geo.pos[ii+1][10][k]; (CE2)
          tordata.uside[ns+1][5][k] =geo.pos[ii+1][11][k]; (CD2)

       and for copyautotogeo - simply rebuild side chain using new chi's (use buildrot routine)
         

      tordata.sidedata[i][0] = geo.iside[ii+1]; //1,2,3,  .. <tordata.sidedata[i][0] = # ofi sidechain units on side i  e.g. if nside =2; then first side is Ca, second is Cb and Cb2
      tordata.sidedata[i][1]= ns;// start number for this chain
  (these are new additions)
      tordata.sidedata[i][2]  = -1 if simply list of C atoms (orig. case)
                              =  rottype  = PHE  (if side chain is phenylanine
      tordata.sidedata[i][3] = 5    = number of atoms in last unit of this chain (after connection atom)
      tordata.sidedata[i][4] = 12  total number of atoms in side chain
  6) Then run dynamics - and find new chi
       Instead of using nextatom to make new sidechain - use the following procedure:
      1) For each unit,  actuallyu rotate the unit by the new chi
          For last unit, rotate all the atoms in the unit
      2) This is done while in runtordyn.  When go back to geo, rebuilt the entire
           side chain using copyautotogeo (buildrot)
  7) I think original findsidechi will work with out any changes

NEW on 23sep98 - Modify so that does not include Cb in force and density
   calculation if geo.cbden = 0 (indicates Cb not in density)
   1) in copygeotoauto - 
NOTE: On 6aug98 - added test so that quits if absolute value of angle
    becomes greater than 2 PI  - 
 NOTE: ON JULY23, 98 - DELETED A LOT OF DEBUGGING STUFF

*/
#include "drawwindow.h"
#include "peptide.h"
#include "refine.h"
#include "trace.h"
#include "tordyn.h"
//#include "menu.h"  // Needed for routine "errormessage"
#include "math.h"
#include "pdbtolev.h"
#include "assignseq.h"
#include "runtordynnew.h"
#include "extendfit.h"

extern struct maiddatafile maidfile;
extern struct tordyn tordata2;
extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version

#if 0
void checkstruc(int **autolist,float ****unit);
void nsidestep(struct geometry& geo,struct tordyn& tordata);
void nextatom(float a1[3],float a2[3],float a3[3],float th,
     float achi,float len,float a[3]);
void nstep1(int de[3][3],int ep[3][3][3],float delt,
    float ****unit,float ****unitv,float *yy,
    float dyy[3],float dr[4][3],float drca[3]);
void makeunitalanine(struct geometry& geo,float ****unit,float ****unitv,
          int ifind,int iend,int iu,float phi,float psi,float delt); 
void makerevunitalanine(struct geometry& geo,float ****unit,float ****unitv,
          int ifind,int iend,int iu,float phi,float psi,float delt); 
void checkstruc(struct geometry& geo,struct tordyn& tordata);
void findangtorque(struct geometry& geo,struct tordyn& tordata,int tcount,
      float hbar[NUNIT][3],float tor[NUNIT][3]);
#endif
#if 0  // Now in runtordynnew.h
void runtordyn(struct griddata& den1,struct geometry& geo,
     struct tordyn& tordata);
void runonesidetordyn(int i,struct griddata& den1,struct geometry& geo,
     struct tordyn& tordata);
#endif
void checkstruc(struct geometry& geo,struct tordyn& tordata);
void findangtorque(struct geometry& geo,struct tordyn& tordata,int tcount,
      float hbar[NUNIT][3],float tor[NUNIT][3]);


extern int verbose; //Defined in trace.c++  if = 1,print some stuff;  = 2 print a lot of stuff



void setupauto(struct geometry& geo,struct tordyn& tordata)
/*
On 7dec98  - Compared new and old torsion dyn routines
Old version:  Calculates velocity for viscous damping by calculation rotation at each step
    The problem with this is that it uses the dq each time which, after long time mayh
    become disconnected from actual atom positions which are determlined using makealanine
New Version:  Keep direct track of velocities with each step.  For some reason, this requre
    different optimal dampinlg  conditions.
  New Version:  best gclin = 1.4gclin (oldversion).
  A vm around 0.3 is best, but seems centered around aabout 2.75 which is now used.
NOTE: DURING THIS TEST - TURNED OFF THE cbden calculation which does not use
    cbden if cb not in density - this seems to make things worse.
  Data used in auto fit routines
NOTE: orig1 and orig2 are no longer used.
  Unit 1 = Ca - Cb
  Unit 2 = Co _ N
  delt=time step
  X[itype][0,1]: for the bond from itype unit to next unit, bond starts at
      atom=0 in itype and ends at atom=1 on next unit
  Returns: natraces  = number of traces of length >aminlen,
      maxaunits  = estimate of the maximum # of units 1 and 2 needed
             =extraf*sum(maintr[i][0]*(delx[0]+delx[1]+dels[2])/3)/)1.5 
           = number of maintr points * corrected for unit length and then
             mult by extrafactor to correct for fact autotrace does not 
             follow maintr.
      maxsunits = estimat of maxinum of side units needed
         maxaunits x 5
NOTE: The big arrays unit, uside and sidedata are allocateda at run time
   and there sie can be increased if needed
      They are allocated with the size tordata.maxaunits and tordata.maxsunits
   The smaller array arrays are preallocated with NUNIT, etc. and ar fixed.
      HOWEVER: tordata.maxaunits must be less than NUNIT1
   Thus, to save space, can initially set  tordata.maxaunits  to small value
     much smaller than NUNIT1 and then increase as needed.
*/
{
    int i;

 tordata.useringang = 0;  //Default - do not allow varing ringangle - if = 1  allow ringangle to move during dynamics
// tordata.useringang = 1;  //if = 1  allow ringangle to move during dynamics
 tordata.ringlim = 10.0*PI/180.0; // all ringang to vary +/- 10 degres
 tordata.delt= 0.002;
    tordata.delt= 0.01;// Works well
 tordata.delt=0.09;// Causes problems
 tordata.delt=0.03;// Works OK but,I think, gives lower density,and does not shorten tiem in batch mode
    tordata.delt= 0.01;// Works well 
//tordata.delt= 0.005;// Works well   
    tordata.useboneconstr = 0;  // default - Do not Use bone constraint on Ca (or N?) in dynamics
    tordata.bonefor = 1.0;  // Force keeping Ca (and N) close to edited bone
    tordata.fixend=0;  // default value, if = 3 add force on end atoms
    tordata.adata[1]=1; /* Indcicate need to  allocate arrays run autoroutine */
    tordata.adata[2]=1; /* Start fit on trace 1 */
    tordata.adata[3] = 30;  // Number of dynamic steps before drawing new peptide
   tordata.vm=0.40; //Used most of summer 98 - Best value for old tordyn (see above)
//tordata.vm= 0.275;  // Slightly better value for new dyn. (see abaove) - used with tordata.glin = 1.4*tordata.gclin
   tordata.gclin=2.0/5.0; /*linear visous damping -best value with old (velocity) tordn (using tordata.vm = 0.3) */
//tordata.gclin= 0.5*tordata.gclin; /* TOO SMALL  linear visous damping */
tordata.gclin=1.6*tordata.gclin;//-best value with old (velocity) tordn (using tordata.vm = 0.3)
    tordata.rfor=sqrt( (2.0/3.0)*tordata.vm*tordata.vm*(tordata.gclin)/(tordata.delt) ); // random force
    tordata.rforside= tordata.rfor/2.0;  //NOTE THIS IS NOT USED NOW random force on a.a. side chains.
    tordata.irepeat = 0; // if =1 then repeat tordyn at lower temp
    tordata.torconst= 1.0/maidfile.denscale; /*Force constant for density force in torsional dyn*/
    tordata.torconstside = tordata.torconst; // force on side 
    tordata.torconstend = tordata.torconst;  // force on  end residue
//    tordata.fuseconst = 1.0/0.001; // Force that keeps two ends togetther when run dynamics in middle Seems OK - reduces error from about 1 to 0.1 Angstrom 
 tordata.fuseconst = 1.0/0.01; //used upto 9jan00-  Force that keeps two ends togetther when run dynamics in middle Seems OK - reduces error from about 1 to 0.1 Angstrom 
// tordata.fuseconst = 10.0; // Force that keeps two ends togetther when run dynamics in middle Seems OK - reduces error from about 1 to 0.1 Angstrom 
///tordata.torconst= 1.5/maidfile.denscale; /*Force constant for density force in torsional dyn*/
#if 0 // for debuggging
tordata.gclin=10*tordata.gclin; /*linear visous damping */
tordata.rfor= 0.0;
#endif
    for(i=0;i<=3;i++)/* assign 4 atoms (Ca, Cb, Co, N dimensionless mass = 1*/
       tordata.amass[i]=1.0;
    tordata.maxaunits = NUNIT1-2;  //Default number - Set to maximum allowed
    tordata.maxsunits = 5*tordata.maxaunits;  // Default - Total # side chain atoms
    tordata.nprelim = 0;  // the number of tordyn steps in tordynbat before start search for optimin fit
#if 0 // this was used up to July 22, 98
    for(i=0;i<NCONSTRAINTS;i++) // Use same angular constaint for all secondary structure
      tordata.angforconst[i] = 1.0;
#endif
#if 0 // this was after  July 22, 98
    for(i=0;i<NCONSTRAINTS;i++) // Use same angular constaint for all secondary structure
      tordata.angforconst[i] = 6.0;
#endif
#if 1 //USED for all good fits up to 8apr99  Use larger value only for helix   after Dec. 8, 98 to 11feb99
    tordata.angforconst[0] = 6.0; // New - use constraint even if 0 constraint
    tordata.angforconst[1] = 12.0;
//    tordata.angforconst[1] = 24.0; // 24apr99  Try using large value for helix
    for(i=2;i<NCONSTRAINTS;i++) // Use smaller alue for the rest
      tordata.angforconst[i] = 6.0; 
#endif
#if 0 // Try large value - Use same angular constaint for all secondary structure
    for(i=0;i<NCONSTRAINTS;i++) 
      tordata.angforconst[i] = 18.0;
#endif
      allocatetordata(tordata);
      tordata.usediamond = 0;  // default, do not use diamond.
//      tordata.usediamond = 1;  // default, do not use diamond.
   if(tordata.usediamond==1){
       printf(" NOTE: USING DIAMOND FORCE \n");
   }
   tordata.nmax = geo.gnmax; // NOTE - this is not used any more??
}

void finddengrad(struct griddata& den1)
/* find gradient of density at all postions
  dengr[i][j][k][ii] is the ii (x,y orz) gradient at i+1/2, j,k
  Note dengr is signed char while den is unsigned. It is possible that difference
   in den is greater then 127, (but unlikely) in which case will get wrong value
   for grad.
*/
{
   int i,j,k;

   for(i=1;i<=den1.amax[0]-den1.amin[0];i++)
     for(j=1;j<=den1.amax[1]-den1.amin[1];j++)
       for(k=1;k<=den1.amax[2]-den1.amin[2];k++){
         den1.dengr[i][j][k][0]=den1.den[i][j][k]-den1.den[i-1][j][k];         
         den1.dengr[i][j][k][1]=den1.den[i][j][k]-den1.den[i][j-1][k];
         den1.dengr[i][j][k][2]=den1.den[i][j][k]-den1.den[i][j][k-1];
       }
}


void findforden(struct griddata& den1,struct geometry& geo,float apos[3],float forden[3])
/*
Modified from old maid routines:
New version (6/25/98)  if point out of range, set geo.quit = 1 - indicates should quit this trace
New version  (10/12/95) if point out of range - do not move back to boundary
    but just pretend it is at boundary when finding force
    That is keep ifix = 0 at all times.
 given atom postion, find force vector due to  density 
 Procedure: e.g. x component:
   grad[i][j][k]=den[i][j][k]-den[i-1][j][k] is actually at position grad[i-1/2][j][k]
 the lower left corner of cube used to interopolate is at round down of (x+1/2). that is
  if x,y,z is (3.7,2.3.4.5) then low. corner is at grad[4][2][4]. When interolate use
   (4.5-x) and (x-3.5) as factors. 
  da[0]=0;da[1]=1;
*/
{
    int i,j,k,ir[3],irh[3];
    float fd[3];

    xtof(den1,apos,fd);
/* New version, interpolates all neighbors*/
    for(k=0;k<=2;k++){
       ir[k]=fd[k]; /*round down*/
       irh[k]=fd[k]+0.5; /*round to next 1/2*/
       if( (ir[k]<1)|| (irh[k]>den1.amax[k]-den1.amin[k]-1) ){
//         if(verbose>=1)
         if(verbose>=2) // 30nov00 - do not print this line
            fprintf(fpout,"ir[%d] = %d irh = %d imax=%d Point too close to boundary in findforden,SETTING geo.quit = 1\n",
                  k,ir[k],irh[k],den1.amax[k]-den1.amin[k]-1);
          geo.quit=1; // indicates should quit this trace
          goto quitfindfor;
       }
       forden[k]=0.0;
    }
    for(i=0;i<2;i++) for(j=0;j<2;j++) for(k=0;k<2;k++) {/*interpolate x,y,z gradient*/
       forden[0]= forden[0]+(int)den1.dengr[irh[0]+geo.da[i]][ir[1]+geo.da[j]][ir[2]+geo.da[k]][0]*
                 fabs( (fd[0]-irh[0]+0.5+geo.al[i])*(fd[1]-ir[1]+geo.al[j])*(fd[2]-ir[2]+geo.al[k]) );
       forden[1]= forden[1]+(int)den1.dengr[ir[0]+geo.da[i]][irh[1]+geo.da[j]][ir[2]+geo.da[k]][1]*
                 fabs( (fd[0]-ir[0]+geo.al[i])*(fd[1]-irh[1]+0.5+geo.al[j])*(fd[2]-ir[2]+geo.al[k]) );
       forden[2]= forden[2]+(int)den1.dengr[ir[0]+geo.da[i]][ir[1]+geo.da[j]][irh[2]+geo.da[k]][2]*
                 fabs( (fd[0]-ir[0]+geo.al[i])*(fd[1]-ir[1]+geo.al[j])*(fd[2]-irh[2]+0.5+geo.al[k]) );
    
    }
quitfindfor:;
}




void allocatetordata(struct tordyn& tordata)
/* maxaunits = max. # of unit1 or unit2 = max # of residues
 NOTE:  This cannot be larger than NUNIT1
 NEW Modified to allow use of a.a. residue
  means tordata.sidedata[i][0 to 3]
  tordata.uside[][12]  // to accomodate tryptophan
*/
{
      fmat4(tordata.unit,2,tordata.maxaunits,4,3);
      fmat4(tordata.unitv,2,tordata.maxaunits,4,3);
      fmat3(tordata.uside,tordata.maxsunits,12,3);
      fmat3(tordata.usidev,tordata.maxsunits,12,3);
      imat2(tordata.sidedata,tordata.maxsunits,5);
}
void deallocatetordata(struct tordyn& tordata)
{
      freefmat4(tordata.unit,2,tordata.maxaunits,4,3);
      freefmat4(tordata.unitv,2,tordata.maxaunits,4,3);
      freefmat3(tordata.uside,tordata.maxsunits,12,3);
      freefmat3(tordata.usidev,tordata.maxsunits,12,3);
      freeimat2(tordata.sidedata,tordata.maxsunits,5);
}

void findtorphipsi(struct tordyn& tordata)
/* Find the phipsi angles for the fit to mainchain itrace,
    In forward direciton:
         psi = rotation between unit1 and unit2 (about Ca - Co bond)
         phi = rotation between unit2 and unit1 (about N - Ca bond)
  unit[0][i][0,1][x,y,z] = x,y,z pos of ith unit1 of Ca and Cb
         if start with unit1, then unit[0][1][2]= N
         if end with unit1 then unit[0][n][2]=Co
     unit[1][i][0,1,2][x,y,z] = x,y,z,pos of ith unit2 of Co  N and O
         if end with unit2, then unit[1][2][n]= Ca
*/
{
   int i;
   int count;

   count=1;
   // find chi =first phi: unit2[i] to unit 1[i]  (unit[1][1][3] = previous Ca)
   i=1;
#if 0
   findchi(tordata.unit[1][i][3],tordata.unit[1][i][1],
                  tordata.unit[0][i][0],tordata.unit[1][i+1][0],&tordata.phipsi[count]);
#endif
#if 1 // gives same result as above since Ca and C in same plane
   findchi(tordata.unit[1][i][0],tordata.unit[1][i][1],
                  tordata.unit[0][i][0],tordata.unit[1][i+1][0],&tordata.phipsi[count]);
   tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
#endif
//fprintf(fpout,"phipsi (degrees) = %6.3f   i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
   for(i=1;i<=tordata.nunit1-1;i++){ // general case
      count = count+1;
      // find chi = psi: unit1[i] to unit 2[i+1]
      findchi(tordata.unit[1][i][1],tordata.unit[0][i][0],
                  tordata.unit[1][i+1][0],tordata.unit[1][i+1][1],&tordata.phipsi[count]);
      tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
//fprintf(fpout,"psi (degrees) = %6.3f    i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
      count = count+1;
      // find chi = phi: unit2[i+1] to unit 1[i+1]
      findchi(tordata.unit[1][i+1][0],tordata.unit[1][i+1][1],
                  tordata.unit[0][i+1][0],tordata.unit[1][i+2][0],&tordata.phipsi[count]);
      tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
//fprintf(fpout,"phi (degrees) = %6.3f    i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
   }
   // find chi =last psi: unit1[i] to unit 2[i+1] 
   count = count+1;
   findchi(tordata.unit[1][i][1],tordata.unit[0][i][0],
                  tordata.unit[1][i+1][0],tordata.unit[1][i+1][1],&tordata.phipsi[count]);
   tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
//fprintf(fpout,"phipsi (degrees) = %6.3f   i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
}

void findrevtorphipsi(struct tordyn& tordata)
/* Find the phipsi angles for the fit to mainchain itrace,
    In backward  direciton:
         psi = rotation between unit2 and unit1 (about Co - Ca bond)
         phi = rotation between unit1 and unit2 (about Ca - N bond)
  unit[0][i][0,1][x,y,z] = x,y,z pos of ith unit1 of Ca and Cb
         if start with unit1, then unit[0][1][2]= N
         if end with unit1 then unit[0][n][2]=Co
     unit[1][i][0,1][x,y,z] = x,y,z,pos of ith unit2 of Co and N
         if end with unit2, then unit[1][2][n]= Ca
*/
{
   int i;
   int count;

   count=1;
   // find chi =first psi in reverse direction: unit2[i] to unit 1[i]  (unit[1][1][3] = previous Ca)
   i=1;
   findchi(tordata.unit[1][i][1],tordata.unit[1][i][0],
                  tordata.unit[0][i][0],tordata.unit[1][i+1][1],&tordata.phipsi[count]);
   tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
//fprintf(fpout,"phipsi (degrees) = %6.3f   i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
   for(i=1;i<=tordata.nunit1-1;i++){ // general case
      count = count+1;
      // find chi = phi: unit1[i] to unit 2[i+1]
      findchi(tordata.unit[1][i][0],tordata.unit[0][i][0],
                  tordata.unit[1][i+1][1],tordata.unit[1][i+1][0],&tordata.phipsi[count]);
      tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
//fprintf(fpout,"phipsi (degrees) = %6.3f   i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
      count = count+1;
      // find chi = psi: unit2[i+1] to unit 1[i+1]
      findchi(tordata.unit[1][i+1][1],tordata.unit[1][i+1][0],
                  tordata.unit[0][i+1][0],tordata.unit[1][i+2][1],&tordata.phipsi[count]);
      tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
//fprintf(fpout,"phipsi (degrees) = %6.3f   i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
   }
   // find chi =last phi: unit1[i] to unit 2[i+1] 
   count = count+1;
   findchi(tordata.unit[1][i][0],tordata.unit[0][i][0],
                  tordata.unit[1][i+1][1],tordata.unit[1][i+1][0],&tordata.phipsi[count]);
   tordata.oldphipsi[count] = tordata.phipsi[count]; // Initially no change in phipsi
//fprintf(fpout,"phipsi (degrees) = %6.3f   i= %d  count = %d\n",tordata.phipsi[count]*180.0/PI,i,count);
}

void findsidechi(int idir,struct tordyn& tordata,struct geometry& geo)
/*
  idir = 0 for N to C direction, 1 for C to N direction
  finds chi[i] angle for side  chain bond between C(uside[i-2][0]) and C(uside[i-1][0]) 
     except for last C which is C(uside[i-2][1]
    NOTE this chi[i] should correspond to ddqs[i] 
  Use routine findchi for old MAID
 NOTE THIS ASSUMED first unit is of type unit2
 ALSO ASSUMES THAT uside[sunit][0] has been set = Cb
*/
{
    int i,k,sunit,lunit,ii;
    float nvec[3];

    for(i=1;i<=tordata.nunit1;i++){ // Set # sidechain atoms = 0
       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 */
          if(idir==0) // N to C direciton
             ii = i;
          else if(idir==1) //C to N direciton
             ii = i+1;
          for(k=0;k<=2;k++) // Prevous N pos.
              nvec[k]=tordata.unit[1][ii][1][k];
//printf("nvec = %5.3f %5.3f %5.3f\n",nvec[0],nvec[1],nvec[2]);
//printf("Ca = %5.3f %5.3f %5.3f\n",tordata.uside[sunit-1][0][0],tordata.uside[sunit-1][0][1],tordata.uside[sunit-1][0][2]);
          if(sunit==lunit) // then only 1 unit consisting of cb and cb2
            findchi(nvec,tordata.uside[sunit-1][0],tordata.uside[sunit][0],tordata.uside[sunit][1],
                       &tordata.chi[sunit-1]);//uside[sunit-1]= = Ca
          else
             findchi(nvec,tordata.uside[sunit-1][0],tordata.uside[sunit][0],tordata.uside[sunit+1][0],
                       &tordata.chi[sunit-1]);
          tordata.oldchi[sunit-1] =tordata.chi[sunit-1];
          if(sunit<lunit){
             if(sunit+1==lunit)
               findchi(tordata.uside[sunit-1][0],tordata.uside[sunit][0],tordata.uside[sunit+1][0],
                        tordata.uside[sunit+1][1],&tordata.chi[sunit]);
             else
               findchi(tordata.uside[sunit-1][0],tordata.uside[sunit][0],tordata.uside[sunit+1][0],
                   tordata.uside[sunit+2][0],&tordata.chi[sunit]);
             tordata.oldchi[sunit] =tordata.chi[sunit];
          }             
          for(k=sunit+2;k<=lunit;k++){
             if(k==lunit)
               findchi(tordata.uside[k-2][0],tordata.uside[k-1][0],tordata.uside[k][0],
                        tordata.uside[k][1],&tordata.chi[k-1]);
             else
               findchi(tordata.uside[k-2][0],tordata.uside[k-1][0],tordata.uside[k][0],
                   tordata.uside[k+1][0],&tordata.chi[k-1]);
             tordata.oldchi[k-1] =tordata.chi[k-1];
          }
#if 1 //INSERT NEW ringangle STUFF HERE:
float uv[3],vv[3],mag1,mag2;
          if( (tordata.sidedata[i][2]==TYR)||(tordata.sidedata[i][2]==HIS)||
                     (tordata.sidedata[i][2]==PHE)||(tordata.sidedata[i][2]==TRP) ){ // if ring compound allow bond angle to  move
             for(k=0;k<=2;k++){
                uv[k] = tordata.uside[lunit-2][0][k]-tordata.uside[lunit-1][0][k]; //lunit = start of last unit - Cb to Ca vector (use i=0 because first atom in unit
                vv[k] = tordata.uside[lunit][0][k]-tordata.uside[lunit-1][0][k]; // Cb to Cb2 vector
             }
             mag1 = sqrt( dot(uv,uv) );
             mag2 = sqrt( dot(vv,vv) );
             tordata.ringang[i]  = acos( dot(uv,vv)/(mag1*mag2) );
             if(tordata.ringang[i]> (geo.phe + tordata.ringlim) )
                      tordata.ringang[i] = geo.phe + tordata.ringlim;
             else if(tordata.ringang[i]< (geo.phe - tordata.ringlim) )
                      tordata.ringang[i] = geo.phe - tordata.ringlim; 
             tordata.oldringang[i] = tordata.ringang[i];
          }// a ring
#endif// END OF NEW INSERTION:     
       } /* if sidedata */
    } /* i */
}                                



void copyautotogeo(struct geometry& geo,struct tordyn& tordata)
/*
New Version - if sidechain is a.a. residue - then rebuild rotomer using buildrot routine
  dir = 0  or  1 for N to C or C to N direction
  NOTE: that geo starts and ends with Ca
    Thus: actually start and end with unit2 (0-C-N) and add Ca as extra atom 
      I think this atom can be anywhere, just added for correct inertia
    In old bonedyn, always assumed start with unit1, so I hope it works
       to start with unit2
  
  Also assume: All unit1s have same number of sidechains atoms - = geo.nside -1
  For unit1 - 0 = Ca, 1 = Cb (Note, if use side chain, then do not use Cb??)
  For unit2   0 = C, 1= N, 2 = 0
  Start with unit2: 1, 2, ... geo.nend[1]-1
    then to unit1:  1,2, ... geo.nend[1]-2
    Also, it seems need to define previoius atom 
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
     
*/
{
   int i,ii,j,k,ns,atomn,rottype;
   float chiang[3],rotpos[NATOMS][3];

   for(i=1;i<=geo.begend[1]-geo.begend[0];i++){
      if(geo.dir ==0){ // N to c direction
         ii = geo.begend[0]+i-1;
         for(k=0;k<=2;k++){
            geo.pos[ii][2][k]=tordata.unit[1][i][0][k]; /* Co pos */
            geo.pos[ii+1][4][k]=tordata.unit[1][i][1][k]; /* N pos of next residue3 */
            geo.pos[ii][3][k]=tordata.unit[1][i][2][k]; /* O pos */
         }             
         if(i<tordata.nunit2){
           for(k=0;k<=2;k++){
             geo.pos[ii+1][1][k]=tordata.unit[0][i][0][k]; /* Ca pos */
             geo.pos[ii+1][5][k]=tordata.unit[0][i][1][k]; /* Cb pos */
           }
         } 
//printf(" cbdist[%d] = %5.3f\n",ii+1,distv(geo.pos[ii+1][5],geo.pos[ii+1][1]) );
      } // dir = 0
      else{ // C to N
         ii =geo.begend[1]-i;
         for(k=0;k<=2;k++){
            geo.pos[ii][2][k]=tordata.unit[1][i][0][k]; /* Co pos */
            geo.pos[ii+1][4][k]=tordata.unit[1][i][1][k]; /* N pos of next residue3 */
            geo.pos[ii][3][k]=tordata.unit[1][i][2][k]; /* O pos */
         }             
         if(i<tordata.nunit2){
           for(k=0;k<=2;k++){
             geo.pos[ii][1][k]=tordata.unit[0][i][0][k]; /* Ca pos */
             geo.pos[ii][5][k]=tordata.unit[0][i][1][k]; /* Cb pos */
           }
         } 
      } // dir = 1
   }  // i

   if(geo.dir ==0){// N to c direction
     for(k=0;k<=2;k++){ // add Ca[1] to first and last unit2
        geo.pos[geo.begend[0]][1][k]=tordata.unit[1][1][3][k];
        geo.pos[geo.begend[1]][1][k]=tordata.unit[1][tordata.nunit2][3][k];
     }
   }
   else{ // C to N
      for(k=0;k<=2;k++){ // add Ca[1] to first and last unit2
         geo.pos[geo.begend[1]][1][k]=tordata.unit[1][1][3][k];
         geo.pos[geo.begend[0]][1][k]=tordata.unit[1][tordata.nunit2][3][k];
      }
   }
  for(i=1;i<=tordata.nunit1;i++){ // Run through each unit1, set side chain date
    if(geo.dir ==0) // N to c direction
         ii = geo.begend[0]+i-1;
    else // C to N
         ii =geo.begend[1]-i-1; // This make ii+1 OK for this direction
    if(tordata.sidedata[i][0]>0){
      ns=tordata.sidedata[i][1];// start number for this chain = Ca  (Note Ca assigned above)
      atomn = 4;  // side chain atoms in geo are numberd 5 (Cb),6,7...
      if(tordata.sidedata[i][2]== -1){ // Then side chain just list of C's, build old way by adding each atom
        for(j=1;j<geo.iside[ii+1];j++){ // geo.nside = # of side chain atoms
           atomn = atomn+1;  //  5,6,7 ...
           ns = ns+1; //next atom (Cb, etc.)
          if(j < geo.iside[ii+1] -1){  //check if last two side atoms   (they go on one unit)
              for(k=0;k<=2;k++)
                 geo.pos[ii+1][atomn][k]=tordata.uside[ns][0][k];
           } 
           else { // last 2 atoms  
              for(k=0;k<=2;k++){
                geo.pos[ii+1][atomn][k]=tordata.uside[ns][0][k] ; 
                geo.pos[ii+1][atomn+1][k]=tordata.uside[ns][1][k]; 
              }
//fprintf(fpout,"i = %d j = %d ns = %d  k=0  tordata.uside[ns][0][k]= %f\n",i,j,ns,tordata.uside[ns][1][0]);
              j = j+1; // Since added two atoms
           }
         }// j
       } // sidedata == -1
       else{  // side chain is a.a residue
//printf("Sidechain is residue\n");
         rottype = tordata.sidedata[i][2];
         for(j=tordata.sidedata[i][1];j<tordata.sidedata[i][1]+tordata.sidedata[i][0]-1;j++)/* start at Cb*/
            chiang[j-tordata.sidedata[i][1]]=tordata.chi[j];
#if 0 //OLD - before 29jan01  - no ringang
         buildrot(geo,geo.pos[ii+1][4],geo.pos[ii+1][1],geo.pos[ii+1][5],
                   rottype,chiang,rotpos);
#endif
#if 1 // New Insertion to allow variable ringangle
         if( (tordata.sidedata[i][2]==TYR)||(tordata.sidedata[i][2]==HIS)||
                     (tordata.sidedata[i][2]==PHE)||(tordata.sidedata[i][2]==TRP) ) // if ring compound allow bond angle to  move
              buildrotring(geo,geo.pos[ii+1][4],geo.pos[ii+1][1],geo.pos[ii+1][5],
                   rottype,chiang,tordata.ringang[i],rotpos);
         else
              buildrot(geo,geo.pos[ii+1][4],geo.pos[ii+1][1],geo.pos[ii+1][5],
                   rottype,chiang,rotpos);
#endif
         copyonerottogeo(rotpos,geo,ii+1,rottype);// copy rotomer side chain to geo
      }
    } // if sidedata>0
  } // i;
}

float aveatomtor(struct griddata& den1,struct geometry& geo,float pos[3])
// Return density for atom postion in tordata co-ord
// Similar to aveatom
{
   int i,j,k,kc,ir[3];
   float fd[3],atomden;

   atomden = 0;
   xtof(den1,pos,fd); /* convert from double word to float sccreen*/
   for(kc=0;kc<=2;kc++)
       ir[kc]=fd[kc]; /*round down to integer*/
   for(k=0;k<=2;k++){ // check if pos in density range
      if( (ir[k]<2)|| (ir[k]>den1.amax[k]-den1.amin[k]-2) ){
         if(verbose>=2)
            fprintf(fpout," ERROR: atom out of range in aveatomtor, returning den = -1000.0 \n");
#if 0 // New 6jun00 - do not set geo.quit=1 if out of range, just return bad density
         geo.quit = 1;
#endif
         return -1000.0;
      }
   }
   for(i=0;i<2;i++) for(j=0;j<2;j++) for(k=0;k<2;k++) {/*interpolate*/
       atomden=atomden+ (int)den1.den[ir[0]+geo.da[i]][ir[1]+geo.da[j]][ir[2]+geo.da[k]]*
                 fabs( (fd[0]-ir[0]+geo.al[i])*(fd[1]-ir[1]+geo.al[j])*(fd[2]-ir[2]+geo.al[k]) );
   }
   return atomden;
}

void avestddevoneside(struct griddata& den1,struct geometry& geo,
      struct tordyn& tordata,int iside,float *aveden,
      int *nbad,int *badset,float *stddev) 
/*
NOTE: NOT USED - IF USED - NEED TO FIX checkatom call conditions
Returns standdev of density and
  Keeps track of bad contact with set
  Does not check for Ala and Pro - just one unit
    and does not check first unit
 find density for side chain = iside
 ns is start, corresponds to Ca, Cb = ns+1, etc.
 if geo.phenylden = 1, then only find density on last unit
*/
{
   int i,j,ns,count;
   float aden[20];

   ns=tordata.sidedata[iside][1];// start number for this chain = Ca
   count = 0;
   *aveden = 0;
   *nbad = 0;
   *badset = 0;
   *stddev = 0;
   geo.quit = 0; // used as marker for atom out range
   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;      
        aden[count]= aveatomtor(den1,geo,tordata.uside[ns+1][i]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto quitrout; 
      }    
   } //Ala and Pro
   else{ // not Ala
     if(geo.phenylden!=1){ // find den of atoms up to last unit
       for(i=ns+1;i<=ns+tordata.sidedata[iside][0]-2;i++){
        count=count+1;
        aden[count]= aveatomtor(den1,geo,tordata.uside[i][0]);  // CB,etc.
        if(geo.quit==1) // means atom out of range
           goto quitrout; 
        if(i>ns+1){ // do not check first unit
           if(checkatom(den1,geo,geo.set,tordata.uside[i][0])== 0) // if atom bad
              *badset = *badset+1;
        }
      }
     } // phenyl
     i=ns+tordata.sidedata[iside][0]-1;
     for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         count=count+1;
         aden[count] = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto quitrout; 
         if(checkatom(den1,geo,geo.set,tordata.uside[i][j])== 0) // if atom bad
              *badset = *badset+1;
     }
   } // not Ala
   for(i=1;i<=count;i++)
      *aveden = *aveden + aden[i];
   *aveden = *aveden/count;
   for(i=1;i<=count;i++)
      *stddev = *stddev + sqr(aden[i]- *aveden);
   *stddev = *stddev/count;
   *stddev = sqrt(*stddev);
   quitrout:;
}  


void avedensetoneside(struct griddata& den1,struct geometry& geo,
      struct tordyn& tordata,int iside,float *aveden,
      int *nbad,int *badset) 
/*
Keeps track of bad contact with set
  Does not check for Ala and Pro - just one unit
    and does not check first unit
 find density for side chain = iside
 ns is start, corresponds to Ca, Cb = ns+1, etc.
 if geo.phenylden = 1, then only find density on last unit
*/
{
   int i,j,ns,count;
   float aden;

   ns=tordata.sidedata[iside][1];// start number for this chain = Ca
   count = 0;
   *aveden = 0;
   *nbad = 0;
   *badset = 0;
   geo.quit = 0; // used as marker for atom out range
   if(tordata.sidedata[iside][0]==0) { // Gly - no unit
      goto quitrout;
   }
   else if(tordata.sidedata[iside][0]==1){ //NOTE- THIS IS NO LONGER CONDITION FOR ALA ANND PRO Ala and Pro - just one unit
      for(i=0;i<=tordata.sidedata[iside][3];i++){
        count=count+1;      
        aden= aveatomtor(den1,geo,tordata.uside[ns+1][i]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto quitrout; 
        *aveden = *aveden+aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
      }    
   } //Ala and Pro
   else{ // not Ala
     if(geo.phenylden!=1){ // find den of atoms up to last unit
       for(i=ns+1;i<=ns+tordata.sidedata[iside][0]-2;i++){
        count=count+1;
        aden= aveatomtor(den1,geo,tordata.uside[i][0]);  // CB,etc.
        if(geo.quit==1) // means atom out of range
           goto quitrout; 
        *aveden = *aveden+aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
        if( i>ns+1){ // do not check first atom
           if(checkatom(den1,geo,geo.set,tordata.uside[i][0])== 0) // if atom bad
              *badset = *badset+1;
        }
      }
     } // phenyl
     i=ns+tordata.sidedata[iside][0]-1;
     for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         count=count+1;
#if 0
printf("i=%d j= %d pos = %5.3f %5.3f %5.3f",
     i,j,tordata.uside[i][j][0],tordata.uside[i][j][1],tordata.uside[i][j][2]);
printf(" done \n");
#endif
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto quitrout; 
         *aveden = *aveden+aden;
         if(aden< geo.minrigid)
           *nbad = *nbad+1;
         if(tordata.sidedata[iside][2]!= PRO){ // do not check proline because second atom is bonded to N
           if(  (tordata.sidedata[iside][0]>2)|| (j>0) ){ // >2  means there are some atoms before lastunit - if =2 then do not use first atom in last unit = Cb
             if(checkatom(den1,geo,geo.set,tordata.uside[i][j])== 0) // if atom bad
               *badset = *badset+1;
           }
         }
     }
   } // not Ala
   *aveden = *aveden/count;
   quitrout:;
}  

void avedenoneside(struct griddata& den1,struct geometry& geo,
      struct tordyn& tordata,int iside,float *aveden,int *nbad) 
// find density for side chain = iside
// ns is start, corresponds to Ca, Cb = ns+1, etc.
// if geo.phenylden = 1, then only find density on last unit
{
   int i,j,ns,count;
   float aden;

   ns=tordata.sidedata[iside][1];// start number for this chain = Ca
   count = 0;
   *aveden = 0;
   *nbad = 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;      
        aden= aveatomtor(den1,geo,tordata.uside[ns+1][i]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto quitrout; 
        *aveden = *aveden+aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
      }    
   } //Ala and Pro
   else{ // not Ala
     if(geo.phenylden!=1){ // find den of atoms up to last unit
       for(i=ns+1;i<=ns+tordata.sidedata[iside][0]-2;i++){
        count=count+1;
        aden= aveatomtor(den1,geo,tordata.uside[i][0]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto quitrout; 
//printf("tor aden = %5.3f\n",aden);
        *aveden = *aveden+aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
      }
     } // phenyl
     i=ns+tordata.sidedata[iside][0]-1;
     for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         count=count+1;
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto quitrout; 
//printf("tor aden = %5.3f\n",aden);
         *aveden = *aveden+aden;
         if(aden< geo.minrigid)
           *nbad = *nbad+1;
     }
   } // not Ala
   *aveden = *aveden/count; 
   quitrout:;
}  
 
void detaildenoneside(struct griddata& den1,struct geometry& geo,
      struct tordyn& tordata,int ibeg,int iend,float *aveden,int *nbad,
      float *firstden,float *lastden,int *totside) 
/* 
NEW 12mar01 - for ILE do not count nbad for i=4 atom = side CH3
  NEW also return geo.res0[1] = nbad in last unit
Identical to avedenoneside except return total aveden,
      den up to last unit and density of last unit
  Useful for ring a.a. and cystein (sulfur atom is last unit)
  totside = # of side atoms
*/
{
   int i,j,ns,count,fcount,lcount,iside;
   float aden;

   count =fcount = lcount= 0;
   *aveden = *lastden = *firstden = 0;
   *nbad = geo.res0[1] =  0;
 for(iside = ibeg;iside<=iend;iside++){
   ns=tordata.sidedata[iside][1];// start number for this chain = Ca
   if(tordata.sidedata[iside][0]==0) { // Gly - no unit
      goto nextiside;
   }
   else if(tordata.sidedata[iside][0]==1){ // NOTE - THIS IS NOT LONGER THE CASE FOR ALA AND PROAla and Pro - just one unit
      for(i=0;i<=tordata.sidedata[iside][3];i++){ // atoms in last unit
        aden= aveatomtor(den1,geo,tordata.uside[ns+1][i]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto nextiside; 
        count=count+1;
        fcount = fcount+1;      
        *aveden = *aveden+aden;
        *firstden = *firstden +aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
      }    
   } //Ala and Pro
   else{ // not Ala
     if(geo.phenylden!=1){ // find den of atoms up to last unit
       for(i=ns+1;i<=ns+tordata.sidedata[iside][0]-2;i++){
        aden= aveatomtor(den1,geo,tordata.uside[i][0]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
        count=count+1;
        fcount = fcount+1;      
        *aveden = *aveden+aden;
        *firstden = *firstden +aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
      }
     } // phenyl
     i=ns+tordata.sidedata[iside][0]-1;
#if 0  // NEW  16jan01 - try commenting out this exception
     if(tordata.sidedata[iside][2]== PRO){ // atoms in last unit - special case, do not find densith for cb2 atom
//printf("SPECIAL CASE - PROLINE\n");
        for(j=0;j<=tordata.sidedata[iside][3];j++){ // tordata.sidedata[iside][3] =2 for Prolineatoms in last unit
         if(j==1){ // do not include this atom
            goto nextatom;
         }
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
         count=count+1;
         lcount = lcount+1;
         *aveden = *aveden+aden;
         *lastden = *lastden + aden;
         if(aden< geo.minrigid){
           geo.res0[1]=geo.res0[1]+1;
           *nbad = *nbad+1;
         }
         nextatom:;
       }//j
     }// if Proline
     else{ // not Proline
#endif
     if(tordata.sidedata[iside][2]== ILE){ // atoms in last unit - special case- for ILE to not count nbad for i=4 atoms
      for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
         count=count+1;
         lcount = lcount+1;
         *aveden = *aveden+aden;
         *lastden = *lastden + aden;
         if( (aden< geo.minrigid)&&(j<tordata.sidedata[iside][3]) ){ // FOR ILE to not count nbad for last atom = side CH3
           geo.res0[1]=geo.res0[1]+1;
           *nbad = *nbad+1;
         }
       }//j
     } // ILE
     else{  // not ILE
      for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
         count=count+1;
         lcount = lcount+1;
         *aveden = *aveden+aden;
         *lastden = *lastden + aden;
         if(aden< geo.minrigid){
           geo.res0[1]=geo.res0[1]+1;
           *nbad = *nbad+1;
         }
       }//j
     } //not  ILE
   } // not Ala
   nextiside:;
 }// iside
 if(count>0)
   *aveden = *aveden/count; 
 else
   *aveden = -100.0;
 if(lcount>0)
     *lastden = *lastden/lcount;
 else
   *lastden = -100.0;
 if(fcount>0)
     *firstden = *firstden/lcount;
 else
   *firstden = -100.0;
 *totside = count;
}  

 
void detaildenonesideset(struct griddata& den1,struct geometry& geo,
      struct tordyn& tordata,int ibeg,int iend,float *aveden,int *nbad,
      float *firstden,float *lastden,int *totside,int icheck) 
/* 
 Identical to detaildenoneside - but checks for bad side contacts.
*/
{
   int i,j,ns,count,fcount,lcount,iside;
   float aden;

   count =fcount = lcount= 0;
   *aveden = *lastden = *firstden = 0;
   *nbad = geo.res0[1] =  0;
 for(iside = ibeg;iside<=iend;iside++){
   ns=tordata.sidedata[iside][1];// start number for this chain = Ca
   if(tordata.sidedata[iside][0]==0) { // Gly - no unit
      goto nextiside;
   }
   else if(tordata.sidedata[iside][0]==1){ // NOTE - THIS IS NOT LONGER THE CASE FOR ALA AND PROAla and Pro - just one unit
      for(i=0;i<=tordata.sidedata[iside][3];i++){ // atoms in last unit
        aden= aveatomtor(den1,geo,tordata.uside[ns+1][i]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto nextiside; 
        count=count+1;
        fcount = fcount+1;      
        *aveden = *aveden+aden;
        *firstden = *firstden +aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
      }    
   } //Ala and Pro
   else{ // not Ala
     if(geo.phenylden!=1){ // find den of atoms up to last unit
       for(i=ns+1;i<=ns+tordata.sidedata[iside][0]-2;i++){
        aden= aveatomtor(den1,geo,tordata.uside[i][0]);  // CB,etc. 
        if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
        count=count+1;
        fcount = fcount+1;      
        *aveden = *aveden+aden;
        *firstden = *firstden +aden;
        if(aden< geo.minrigid)
          *nbad = *nbad+1;
#if 1 //18jan01  New check atoms for bad contact
        if( (iside>=icheck)&&(i>ns+1) ){ // do not check first atom
           if(checkatom(den1,geo,geo.set,tordata.uside[i][0])== 0) // if atom bad
              *nbad = *nbad+10;// make very bad if bad contact
        }
#endif
      }
     } // phenyl
     i=ns+tordata.sidedata[iside][0]-1;
    if(tordata.sidedata[iside][2]== PRO){ //  NEW  9feb01 - Used slighty different version of Proline special caseatoms in last unit - special case, allower lower bad density for cb2 atom
//printf("SPECIAL CASE - PROLINE\n");
       for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
         count=count+1;
         lcount = lcount+1;
         *aveden = *aveden+aden;
         *lastden = *lastden + aden;
         if(j==1){ // allower lower density for this atom
           if(aden< 0.0){
              geo.res0[1]=geo.res0[1]+1;
             *nbad = *nbad+1;
           }
         }
         else{
           if(aden< geo.minrigid){
              geo.res0[1]=geo.res0[1]+1;
             *nbad = *nbad+1;
           }
        }
       }//j
     }// if Proline
     else if(tordata.sidedata[iside][2]== ILE){ //NOTE: this exception is not used since nbad for ILE not used atoms in last unit - special case- for ILE to not count nbad for i=4 atoms
      for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
         count=count+1;
         lcount = lcount+1;
         *aveden = *aveden+aden;
         *lastden = *lastden + aden;
         if( (aden< geo.minrigid)&&(j<tordata.sidedata[iside][3]) ){ // FOR ILE to not count nbad for last atom = side CH3
           geo.res0[1]=geo.res0[1]+1;
           *nbad = *nbad+1;
         }
       }//j
     } // ILE
     else{  // not ILE or PRO
      for(j=0;j<=tordata.sidedata[iside][3];j++){ // atoms in last unit
         aden = aveatomtor(den1,geo,tordata.uside[i][j]);
         if(geo.quit==1) // means atom out of range
           goto nextiside; 
//printf("tor aden = %5.3f\n",aden);
         count=count+1;
         lcount = lcount+1;
         *aveden = *aveden+aden;
         *lastden = *lastden + aden;
         if(aden< geo.minrigid){
           geo.res0[1]=geo.res0[1]+1;
           *nbad = *nbad+1;
         }
#if 1 //18jan00  New check atoms for bad contact
         if( (iside>=icheck)&&(tordata.sidedata[iside][2]!= PRO) ){ // do not check proline because second atom is bonded to N
           if(  (tordata.sidedata[iside][0]>2)|| (j>0) ){ // >2  means there are some atoms before lastunit - if =2 then do not use first atom in last unit = Cb
             if(checkatom(den1,geo,geo.set,tordata.uside[i][j])== 0) // if atom bad
               *nbad = *nbad+10;// make very bad if bad contact
           }
         }
#endif 
      }//j
     } // not ILE
   } // not Ala
   nextiside:;
 }// iside
 if(count>0)
   *aveden = *aveden/count; 
 else
   *aveden = -100.0;
 if(lcount>0)
     *lastden = *lastden/lcount;
 else
   *lastden = -100.0;
 if(fcount>0)
     *firstden = *firstden/lcount;
 else
   *firstden = -100.0;
 *totside = count;
}  


void copyonesidetogeo(struct geometry& geo,struct tordyn& tordata)
/*
    Copies the side resdiue from auto to geo   
  Just the reverse of this part of copygeotoauto
*/
{
   int i,ii,jj,k,ns,rottype;
   
   i=1;
   ns=tordata.sidedata[i][1];// start number for this chain
   ii = geo.begend[0]+i-1;
   for(k=0;k<=2;k++)
           geo.pos[ii+1][1][k]=tordata.uside[ns][0][k]; // first atom is Ca
//printf("Copy residue\n");
      ns = ns+1;
      rottype = -geo.iside[ii+1];
      switch(rottype){
        case(ALA): /* just use the above ideal CB */
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          break;
        case(ARG): //branched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k] ;// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k] ;//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][7][k]=tordata.uside[ns][0][k] ;//
            geo.pos[ii+1][8][k]=tordata.uside[ns][1][k] ;//
            geo.pos[ii+1][9][k]=tordata.uside[ns][2][k] ;//
            geo.pos[ii+1][10][k]=tordata.uside[ns][3][k] ;
            geo.pos[ii+1][11][k]= geo.pos[ii+1][9][k];//
            geo.pos[ii+1][12][k]=tordata.uside[ns][4][k] ;//
          } // k
          break;
        case(ASN): //branched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][1][k];//
             geo.pos[ii+1][8][k]= geo.pos[ii+1][6][k];
            geo.pos[ii+1][9][k]=tordata.uside[ns][2][k];
          } // k
          break;
        case(ASP): //branched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][1][k];//
             geo.pos[ii+1][8][k]= geo.pos[ii+1][6][k];
            geo.pos[ii+1][9][k]=tordata.uside[ns][2][k];
          } // k
          break;
        case(CYS): //branched
          for(k=0;k<=2;k++){
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];//Cb
            geo.pos[ii+1][6][k]=tordata.uside[ns][1][k];//
          } // k
          break;
        case(GLU): //branched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][7][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][8][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][9][k]= geo.pos[ii+1][7][k];
            geo.pos[ii+1][10][k]=tordata.uside[ns][2][k];
          } // k
          break;
        case(AA4): //Artificial straight 4 atoms res- 
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][7][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][8][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][9][k]= geo.pos[ii+1][7][k];
//            geo.pos[ii+1][10][k]=tordata.uside[ns][2][k];
          } // k
          break;
         case(GLN): //branched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][7][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][8][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][9][k]= geo.pos[ii+1][7][k];
            geo.pos[ii+1][10][k]=tordata.uside[ns][2][k];
          } // k
          break;
        case(GLY): 
           break;
        case(HIS): // unbranched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
            for(jj=1;jj<=tordata.sidedata[i][3];jj++)
               geo.pos[ii+1][jj+6][k]=tordata.uside[ns][jj][k];// (CD1)
            geo.pos[ii+1][11][k]= geo.pos[ii+1][6][k];
          } // k
         break;
#if 0 // Old version - ERROR- cannot just omit CH3, because does not know where to put it after dynamics
        case(ILE): //SLIGHT ERROR: DO NOT INCLUDE SIDE CH3 IN dynamic force calculation (ILE does not fit assume rule that only last group has more than one atom
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][1][k];//
          } // k
          break;
#endif
#if 1 // NEW  version -treat all terminal atoms as one rigid unit - neglects one degree of freedom
        case(ILE): //SLIGHT ERROR: DO NOT INCLUDE SIDE CH3 IN dynamic force calculation (ILE does not fit assume rule that only last group has more than one atom
          for(k=0;k<=2;k++){
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
            geo.pos[ii+1][6][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][2][k];// This is far end atom
            geo.pos[ii+1][8][k]=tordata.uside[ns][3][k];//
          } // k
          break;
#endif
        case(LEU): //Same as ASP
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][1][k];//
//            geo.pos[ii+1][8][k]= geo.pos[ii+1][5][k];
            geo.pos[ii+1][8][k]= geo.pos[ii+1][6][k];
            geo.pos[ii+1][9][k]=tordata.uside[ns][2][k];
          } // k
          break;
        case(AA3): //Artificial straifght 3 residue 
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][8][k]= geo.pos[ii+1][6][k];
//            geo.pos[ii+1][9][k]=tordata.uside[ns][2][k];
          } // k
          break;
        case(LYS): //unbranched
         for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][7][k]=tordata.uside[ns][0][k];//
            geo.pos[ii+1][8][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][9][k]=tordata.uside[ns][2][k];//
          } // k
          break;
        case(MET): //unbranched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];// (Cg)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][7][k]=tordata.uside[ns][0][k];//SD
            geo.pos[ii+1][8][k]=tordata.uside[ns][1][k];//CE
          } // k
         break;
        case(PHE): // unbranched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
            for(jj=1;jj<=tordata.sidedata[i][3];jj++)
               geo.pos[ii+1][jj+6][k]=tordata.uside[ns][jj][k];// (CD1)
            geo.pos[ii+1][12][k]= geo.pos[ii+1][6][k];
          } // k
          break;
        case(PRO): 
          for(k=0;k<=2;k++){
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
            geo.pos[ii+1][6][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][2][k];//
          } // k
          break;
        case(SER): 
          for(k=0;k<=2;k++){
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
            geo.pos[ii+1][6][k]=tordata.uside[ns][1][k];//
          } // k
          break;
       case(THR): 
          for(k=0;k<=2;k++){
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
            geo.pos[ii+1][6][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][2][k];//NORE CB not countd
          } // k
          break;
        case(TRP): // unbranched
         for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
            for(jj=1;jj<=8;jj++)
               geo.pos[ii+1][jj+6][k]=tordata.uside[ns][jj][k];// (CD1)
            geo.pos[ii+1][15][k]= geo.pos[ii+1][6][k];
            geo.pos[ii+1][16][k]= geo.pos[ii+1][14][k];
            geo.pos[ii+1][17][k]= geo.pos[ii+1][9][k];
          } // k
          break;
        case(TYR): // unbranched
          for(k=0;k<=2;k++)
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            geo.pos[ii+1][6][k]=tordata.uside[ns][0][k];//(CG)
            for(jj=1;jj<=4;jj++)
               geo.pos[ii+1][jj+6][k]=tordata.uside[ns][jj][k];// (CD1)
            geo.pos[ii+1][11][k]= geo.pos[ii+1][6][k];
            geo.pos[ii+1][12][k]=tordata.uside[ns][5][k];// (CD1)
            geo.pos[ii+1][13][k]=tordata.uside[ns][6][k];// (CD1)
            geo.pos[ii+1][14][k]= geo.pos[ii+1][9][k];
          } // k
          break;
       case(VAL): // Same as THR
         for(k=0;k<=2;k++){
            geo.pos[ii+1][5][k]=tordata.uside[ns][0][k];// (Cb)
            geo.pos[ii+1][6][k]=tordata.uside[ns][1][k];//
            geo.pos[ii+1][7][k]=tordata.uside[ns][2][k];//NORE CB not countd
          } // k
          break;
        default:
          printf("No sidechain data for rottype = %d\n",rottype);
          break;
      }
}


void ncopygeotoauto(struct geometry& geo,struct tordyn& tordata)
/*
New - allow side chains to be amino acid resides
 idicated by geo.iside = -rottype
 Use brute force procedure - for individual residues
 
  dir = 0:  N to C direction;  = 1  C to N
  NOTE:  geo starts and ends with Ca
    Thus: actually start and end with unit2 (0-C-N) and add Ca as extra atom 
      I think this atom can be anywhere, just added for correct inertia
    In old bonedyn, always assumed start with unit1, so I hope it works
       to start with unit2
  For unit1 - 0 = Ca, 1 = Cb (Note, if use side chain, then do not use Cb??)
  For unit2   0 = C, 1= N, 2 = 0
  Start with unit2: 1, 2, ... geo.nend[1]-1
    then to unit1:  1,2, ... geo.nend[1]-2
    Also, it seems need to define previoius atom 
  Also assume: All unit1s have same number of sidechains atoms - = geo.nside -1
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
 uside[i][0,1][x,y,z] = xyz postion of C (0) of side chain i. 
    at End of side chain - have C - C (2 atoms) correspond to 0, and 1. 
 sidedata[i]0]= # of side chains on ith unit1.
 sidedata[i][1]= the # of uside that starts on ith unit1.
 sidedata[0][0]= tot. # of usides that have been currently assigned
 ttype[i] = constraint type for bond preceding unit[i] (of type1 or 2)
Details for sidechains:
  for each side chain (iside) with more than just Cb atom:
    1)i= sidedata[0][0]=sidedata[0][0]+1; // coounts total number of side chians
      sidedata[iside][1]= i;// start number for this chain
      uside[i][0][k] = ca[k] // first atom is Ca
      uside[i+1][0][k] = Cb  // second is Cb 
      sidedta[iside][0] = 2; // # of side chain atoms on this chain
    2) for each additional side atom
      i = sidedata[0][0] = sidedata[0][0]+1; // increase total number by 1
      uside[i][0][k] = atom pos. // assign postion
      sidedata[iside][0] = sidedata[iside][0]+1; // increase numbe on this chain by 1
    3) Then combine last two side chains into 1 unit with two atoms
      i = sidedata[0][0] = sidedata[0][0]-1; // reduce total # of side atoms by 1
      sidedata[iside][0] = sidedata[iside][0]-1;// reduce # of side atoms on this chain by 1
       uside[i][1][k] = uside[i+1][0][k];  // set position of secnd atom on unit
*/
{
   int i,ii,j,jj,k,ns,atomn,count,rottype,tcount;
   if(geo.dir==0){  // N to c direction
      tordata.cbond[0][0]=0;tordata.cbond[0][1]=0; /* bond from unit1 to unit 2 is from Ca to C */
      tordata.cbond[1][0]=1;tordata.cbond[1][1]=0; /* bond from unit2 to unit1 is from N to Ca*/
   }
   else if(geo.dir==1){  // C to N direction
      tordata.cbond[0][0]=0;tordata.cbond[0][1]=1; /* bond from unit1 to unit 2 is from Ca to N */
      tordata.cbond[1][0]=0;tordata.cbond[1][1]=0; /* bond from unit2 to unit1 is from C to Ca*/
   }
   if(geo.begend[1]-geo.begend[0]>NUNIT1-5){
//       errormessage(toplevel,"ERROR: More then NUNIT1 residues, truncating auto");
       fprintf(fpout,"ERROR: More then NUNIT1 residues, Increase NUNIT1 in tordyn.h and recompile\n");  // Cannot use errormessage becuase this is called in batch mode
       geo.quit = 1;
       goto quitrout;
   }
   else{ // Use complete lenght     
      tordata.nunit1=geo.begend[1]-geo.begend[0]-1;
      tordata.nunit2 = geo.begend[1]-geo.begend[0];
   }
   tcount = tordata.nunit1+tordata.nunit2; 
   for(i=1;i<tcount;i++) //NEW 18feb99 - set tordata.angsign[i] =0 so will find correct sign
      tordata.angsign[i]=0; 
   count = 0; // counts unit's (1 and 2)
   for(i=1;i<=geo.begend[1]-geo.begend[0];i++){
      if(geo.dir ==0){ // N to c direction
         ii = geo.begend[0]+i-1;
         tordata.cis[i] = geo.cis[ii];  // set cis corresponding to Ca
#if 0
if(tordata.cis[i]==1)
  printf(" For i = %d, using cis tordata\n",i);
#endif
//fprintf(fpout,"geo.rtype[%d] = %d\n",ii,geo.rtype[ii]);
         count = count+1;
//         tordata.ttype[count]= geo.rtype[ii];  //THIS IS OLD VERSION (1SEP98) AND (I think) is wrong constaint for bond preceding type 2 unit (NOTE: first unit2 has no constraint - give arbitrary value
         tordata.ttype[count]= geo.rtype[ii+1];  // NEW constaint for bond preceding type 1 (=phi) unit (NOTE: first unit2 has no constraint - give arbitrary value
         count = count+1;
         tordata.ttype[count]= geo.rtype[ii+1];  // NEW constaint for bond preceding type 2 (=psi) unit
         for(k=0;k<=2;k++){
           tordata.unit[1][i][0][k]=geo.pos[ii][2][k]; /* Co pos */
           tordata.unit[1][i][1][k]=geo.pos[ii+1][4][k]; /* N pos of next residue3 */
           tordata.unit[1][i][2][k]=geo.pos[ii][3][k]; /* O pos */
         }             
         if(i<geo.begend[1]-geo.begend[0]){
           for(k=0;k<=2;k++){
             tordata.unit[0][i][0][k]=geo.pos[ii+1][1][k]; /* Ca pos */
             tordata.unit[0][i][1][k]=geo.pos[ii+1][5][k]; /* Cb pos */
           }
         } 
      } // dir = 0
      else {// C to N
         ii =geo.begend[1]-i;
         tordata.cis[i] = geo.cis[ii];  // set cis corresponding to Ca
#if 0
if(tordata.cis[i]==1)
  printf(" For i = %d, using cis tordata\n",i);
#endif
         count = count+1;
         tordata.ttype[count]= geo.rtype[ii];  // constaint for bond preceding type 2 unit
         count = count+1;
         tordata.ttype[count]= geo.rtype[ii];  // constaint for bond preceding type 1unit
         for(k=0;k<=2;k++){
           tordata.unit[1][i][0][k]=geo.pos[ii][2][k]; /* Co pos */
           tordata.unit[1][i][2][k]=geo.pos[ii][3][k]; /* O pos */
           tordata.unit[1][i][1][k]=geo.pos[ii+1][4][k]; /* N pos of next residue3 */
         }             
         if(i<geo.begend[1]-geo.begend[0]){
           for(k=0;k<=2;k++){
             tordata.unit[0][i][0][k]=geo.pos[ii][1][k]; /* Ca pos */
             tordata.unit[0][i][1][k]=geo.pos[ii][5][k]; /* Cb pos */
           }
         } 
      } // dir = 1
   }  // i
   if( (geo.fixend==3)&&(geo.dir==0) )// special case where cis proline is the pick point for local tor. dyn
       tordata.cis[geo.begend[1]-geo.begend[0]+1] = geo.cis[ geo.begend[1] ];
   else
       tordata.cis[geo.begend[1]-geo.begend[0]+1] =0;  // Need to define for imax+1
   if(geo.dir ==0){// N to c direction
     for(k=0;k<=2;k++){ // add Ca[1] to first and last unit2
        tordata.unit[1][1][3][k]=geo.pos[geo.begend[0]][1][k];
        tordata.unit[1][tordata.nunit2][3][k]=geo.pos[geo.begend[1]][1][k];
     }
   }
   else{ // C to N
      for(k=0;k<=2;k++){ // add Ca[1] to first and last unit2
         tordata.unit[1][1][3][k]=geo.pos[geo.begend[1]][1][k];
         tordata.unit[1][tordata.nunit2][3][k]=geo.pos[geo.begend[0]][1][k];
      }
   }
   if(geo.dir ==0){
      findtorphipsi(tordata);
      findangsign(geo,tordata); // finds sign of torsional force - if needed for phi/psi constraint
   }
   else if(geo.dir ==1){
      findrevtorphipsi(tordata);
      findangsign(geo,tordata); // finds sign of torsional force - if needed for phi/psi constraint
   }
  ns=tordata.sidedata[0][0] = 0; // counts total number of side chain atoms
  for(i=1;i<=tordata.nunit1;i++){ // Run through each unit1, set side chain date
   if(geo.dir ==0) // N to c direction
         ii = geo.begend[0]+i-1;
   else // C to N
         ii =geo.begend[1]-i-1; // This make ii+1 OK for this direction
   if(geo.iside[ii+1] ==1){ // if only one side chain atom (= Cb, then no atoms in uside
        tordata.sidedata[i][0] = 0;
        tordata.sidedata[i][2] = -1; // Indicate geo.iside >0
        tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
        goto noside; 
   } 
//fprintf(fpout,"ii=%d geo.iside[ii+1]=%d\n",ii,geo.iside[ii+1]);
   ns=ns+1; 
   tordata.sidedata[i][1]= ns;// start number for this chain
   for(k=0;k<=2;k++)
           tordata.uside[ns][0][k] = geo.pos[ii+1][1][k]; // first atom is Ca
   atomn = 4;  // side chain atoms in geo are numberd 5 (Cb),6,7...
   if(geo.iside[ii+1]>0){ // Then standard old proceudre - just sequence of C
        tordata.sidedata[i][0] = geo.iside[ii+1]; //# of side chain atoms e.g. if nside =2; then first side is Ca, second is Cb and Cb2
        tordata.sidedata[i][2] = -1; // Indicate geo.iside >0
        tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
        for(j=1;j<geo.iside[ii+1];j++){ // geo.nside = # of side chain atoms
           atomn = atomn+1;  //  5,6,7 ...
           ns = ns+1; //next atom (Cb, etc.)
           if(j < geo.iside[ii+1] -1){  //check if last two side atoms   (they go on one unit)
              for(k=0;k<=2;k++)
                tordata.uside[ns][0][k] =geo.pos[ii+1][atomn][k];
           } 
           else { // last 2 atoms  
              for(k=0;k<=2;k++){
                tordata.uside[ns][0][k] =geo.pos[ii+1][atomn][k]; 
                tordata.uside[ns][1][k] =geo.pos[ii+1][atomn+1][k]; 
              }
              j = j+1; // Since added two atoms
           } // else
        }// j
    } // iside>=0
    else {// geo.iside <=0
//printf("Copy residue\n");
      ns = ns+1;
      rottype = -geo.iside[ii+1];
      tordata.sidedata[i][2] = rottype; // indicate a.a. residue
      switch(rottype){
        case(ALA): /* just use the above ideal CB */
          tordata.sidedata[i][0] = 2; //NOTE: MUST MAKE THIS 2 FOR runtordyn to work OK, # of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 0;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          tordata.sidedata[i][4] = 5;// Total number of side chain atoms
          break;
        case(ARG): //branched
          tordata.sidedata[i][0] = 4; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 4;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][8][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][9][k];//
            tordata.uside[ns][3][k] =geo.pos[ii+1][10][k];
            tordata.uside[ns][4][k] =geo.pos[ii+1][12][k];//
          } // k
          tordata.sidedata[i][4] = 12;// Total number of side chain atoms
          break;
        case(ASN): //branched
          tordata.sidedata[i][0] = 3; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][9][k];
          } // k
          tordata.sidedata[i][4] = 9;// Total number of side chain atoms
          break;
        case(ASP): //branched
          tordata.sidedata[i][0] = 3; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][9][k];
          } // k
          tordata.sidedata[i][4] = 9;// Total number of side chain atoms
          break;
        case(CYS): //branched
          tordata.sidedata[i][0] = 2; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];//Cb
            tordata.uside[ns][1][k] =geo.pos[ii+1][6][k];//
          } // k
          tordata.sidedata[i][4] = 6;// Total number of side chain atoms
          break;
        case(GLU): //branched
          tordata.sidedata[i][0] = 4; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][8][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][10][k];
          } // k
          tordata.sidedata[i][4] = 10;// Total number of side chain atoms
          break;
        case(AA4): //artificial 4 atom strgint
          tordata.sidedata[i][0] = 4; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][8][k];//
//            tordata.uside[ns][2][k] =geo.pos[ii+1][10][k];
          } // k
          tordata.sidedata[i][4] = 9;// Total number of side chain atoms
          break;
        case(GLN): //branched
          tordata.sidedata[i][0] = 4; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][8][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][10][k];
          } // k
          tordata.sidedata[i][4] = 10;// Total number of side chain atoms
          break;
        case(GLY): 
          tordata.sidedata[i][0] = 0; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 0;//    = number of atoms in last unit of this chain (after connection atom)
          tordata.sidedata[i][4] = 4;// Total number of side chain atoms
          break;
        case(HIS): // unbranched
          tordata.sidedata[i][0] = 3; //# of side chain  units +1e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 4;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
            for(jj=1;jj<=tordata.sidedata[i][3];jj++)
               tordata.uside[ns][jj][k] =geo.pos[ii+1][jj+6][k];// (CD1)
          } // k
           tordata.sidedata[i][4] = 10;// Total number of side chain atoms
         break;
#if 0
        case(ILE): //SLIGHT ERROR: DO NOT INCLUDE SIDE CH3 IN dynamic force calculation (ILE does not fit assume rule that only last group has more than one atom
          tordata.sidedata[i][0] = 3; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][7][k];//
          } // k
          tordata.sidedata[i][4] = 7;// Total number of side chain atoms
          break;
#endif
#if 1 // NEW  version -treat all terminal atoms as one rigid unit - neglects one degree of freedom
        case(ILE): //SLIGHT ERROR: DO NOT INCLUDE SIDE CH3 IN dynamic force calculation (ILE does not fit assume rule that only last group has more than one atom
          tordata.sidedata[i][0] = 2; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 3;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          for(k=0;k<=2;k++){
            tordata.uside[ns][1][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][3][k] =geo.pos[ii+1][8][k];//
          } // k
#if 0
printf("dist cb to 6 = %5.3f\n",distv(geo.pos[ii+1][5],geo.pos[ii+1][6]));
printf("dist cb to 7 = %5.3f\n",distv(geo.pos[ii+1][5],geo.pos[ii+1][7]));
printf("dist cb to 8 = %5.3f\n",distv(geo.pos[ii+1][5],geo.pos[ii+1][8]));
#endif
          tordata.sidedata[i][4] = 8;// Total number of side chain atoms
          break;
#endif
        case(LEU): //Same as ASP
          tordata.sidedata[i][0] = 3; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][9][k];
          } // k
          tordata.sidedata[i][4] = 9;// Total number of side chain atoms
          break;
        case(AA3): //Artificial 3 atom straignt
          tordata.sidedata[i][0] = 3; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][7][k];//
//            tordata.uside[ns][2][k] =geo.pos[ii+1][9][k];
          } // k
          tordata.sidedata[i][4] = 8;// Total number of side chain atoms
          break;
        case(LYS): //unbranched
          tordata.sidedata[i][0] = 4; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][7][k];//
            tordata.uside[ns][1][k] =geo.pos[ii+1][8][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][9][k];//
          } // k
          tordata.sidedata[i][4] = 9;// Total number of side chain atoms
          break;
        case(MET): //unbranched
          tordata.sidedata[i][0] = 4; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];// (Cg)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][7][k];//SD
            tordata.uside[ns][1][k] =geo.pos[ii+1][8][k];//CE
//            tordata.uside[ns][2][k] =geo.pos[ii+1][8][k];//
          } // k
          tordata.sidedata[i][4] = 8;// Total number of side chain atoms
          break;
        case(PHE): // unbranched
          tordata.sidedata[i][0] = 3; //# of side chain  units +1e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 5;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
            for(jj=1;jj<=tordata.sidedata[i][3];jj++)
               tordata.uside[ns][jj][k] =geo.pos[ii+1][jj+6][k];// (CD1)
          } // k
          tordata.sidedata[i][4] = 11;// Total number of side chain atoms
          break;
        case(PRO): 
          tordata.sidedata[i][0] = 2; //NOTE: MUST MAKE THIS 2 FOR runtordyn to work OK, # of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
            tordata.uside[ns][1][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][7][k];//
          } // k
          tordata.sidedata[i][4] = 7;// Total number of side chain atoms
          break;
        case(SER): 
          tordata.sidedata[i][0] = 2; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 1;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
            tordata.uside[ns][1][k] =geo.pos[ii+1][6][k];//
          } // k
          tordata.sidedata[i][4] = 6;// Total number of side chain atoms
          break;
       case(THR): 
          tordata.sidedata[i][0] = 2; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
            tordata.uside[ns][1][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][7][k];//NORE CB not countd
          } // k
          tordata.sidedata[i][4] = 7;// Total number of side chain atoms
          break;
        case(TRP): // unbranched
          tordata.sidedata[i][0] = 3; //# of side chain  units +1e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 8;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
            for(jj=1;jj<=8;jj++)
               tordata.uside[ns][jj][k] =geo.pos[ii+1][jj+6][k];// (CD1)
          } // k
          tordata.sidedata[i][4] = 14;// Total number of side chain atoms
          break;
        case(TYR): // unbranched
          tordata.sidedata[i][0] = 3; //# of side chain  units +1e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 6;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++)
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
          ns=ns+1;
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][6][k];//(CG)
            for(jj=1;jj<=4;jj++)
               tordata.uside[ns][jj][k] =geo.pos[ii+1][jj+6][k];// (CD1)
            tordata.uside[ns][5][k] =geo.pos[ii+1][12][k];// (CD1)
            tordata.uside[ns][6][k] =geo.pos[ii+1][13][k];// (CD1)
          } // k
          tordata.sidedata[i][4] = 13;// Total number of side chain atoms
          break;
       case(VAL): // Same as THR
          tordata.sidedata[i][0] = 2; //# of side chain  units +1 e.g. if nside =2; then first side is Ca, second is Cb and Cb2
          tordata.sidedata[i][3] = 2;//    = number of atoms in last unit of this chain (after connection atom)
          for(k=0;k<=2;k++){
            tordata.uside[ns][0][k] =geo.pos[ii+1][5][k];// (Cb)
            tordata.uside[ns][1][k] =geo.pos[ii+1][6][k];//
            tordata.uside[ns][2][k] =geo.pos[ii+1][7][k];//NORE CB not countd
          } // k
          tordata.sidedata[i][4] = 7;// Total number of side chain atoms
          break;
        default:
          printf("No sidechain data for rottype = %d\n",rottype);
          break;
      }
    } // else
    noside:;
  } //i
  tordata.sidedata[0][0]= ns; // total # side units
  findsidechi(geo.dir,tordata,geo);
  quitrout:;
}


void copygeotoauto(struct geometry& geo,struct tordyn& tordata)
/*
NEW 18feb99 - set angsign = 0 for all angles - this forces findangsign to
    directly calculate value of sign
  dir = 0:  N to C direction;  = 1  C to N
  NOTE:  geo starts and ends with Ca
    Thus: actually start and end with unit2 (0-C-N) and add Ca as extra atom 
      I think this atom can be anywhere, just added for correct inertia
    In old bonedyn, always assumed start with unit1, so I hope it works
       to start with unit2
  For unit1 - 0 = Ca, 1 = Cb (Note, if use side chain, then do not use Cb??)
  For unit2   0 = C, 1= N, 2 = 0
  Start with unit2: 1, 2, ... geo.nend[1]-1
    then to unit1:  1,2, ... geo.nend[1]-2
    Also, it seems need to define previoius atom 
  Also assume: All unit1s have same number of sidechains atoms - = geo.nside -1
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
 uside[i][0,1][x,y,z] = xyz postion of C (0) of side chain i. 
    at End of side chain - have C - C (2 atoms) correspond to 0, and 1. 
 sidedata[i]0]= # of side chains on ith unit1.
 sidedata[i][1]= the # of uside that starts on ith unit1.
 sidedata[0][0]= tot. # of usides that have been currently assigned
 ttype[i] = constraint type for bond preceding unit[i] (of type1 or 2)
Details for sidechains:
  for each side chain (iside) with more than just Cb atom:
    1)i= sidedata[0][0]=sidedata[0][0]+1; // coounts total number of side chians
      sidedata[iside][1]= i;// start number for this chain
      uside[i][0][k] = ca[k] // first atom is Ca
      uside[i+1][0][k] = Cb  // second is Cb 
      sidedta[iside][0] = 2; // # of side chain atoms on this chain
    2) for each additional side atom
      i = sidedata[0][0] = sidedata[0][0]+1; // increase total number by 1
      uside[i][0][k] = atom pos. // assign postion
      sidedata[iside][0] = sidedata[iside][0]+1; // increase numbe on this chain by 1
    3) Then combine last two side chains into 1 unit with two atoms
      i = sidedata[0][0] = sidedata[0][0]-1; // reduce total # of side atoms by 1
      sidedata[iside][0] = sidedata[iside][0]-1;// reduce # of side atoms on this chain by 1
       uside[i][1][k] = uside[i+1][0][k];  // set position of secnd atom on unit
*/
{
   int i,ii,j,k,ns,atomn,count,tcount;

   if(geo.dir==0){  // N to c direction
      tordata.cbond[0][0]=0;tordata.cbond[0][1]=0; /* bond from unit1 to unit 2 is from Ca to C */
      tordata.cbond[1][0]=1;tordata.cbond[1][1]=0; /* bond from unit2 to unit1 is from N to Ca*/
   }
   else if(geo.dir==1){  // C to N direction
      tordata.cbond[0][0]=0;tordata.cbond[0][1]=1; /* bond from unit1 to unit 2 is from Ca to N */
      tordata.cbond[1][0]=0;tordata.cbond[1][1]=0; /* bond from unit2 to unit1 is from C to Ca*/
   }
   if(geo.begend[1]-geo.begend[0]>NUNIT1-5){
//       errormessage(toplevel,"ERROR: More then NUNIT1 residues, truncating auto");
       fprintf(fpout,"ERROR: More then NUNIT1 residues, truncating auto\n");  // Cannot use errormessage becuase this is called in batch mode
       tordata.nunit1=NUNIT1-2; // Truncate length of auto
       tordata.nunit2 = NUNIT1-2;
   }
   else{ // Use complete lenght     
      tordata.nunit1=geo.begend[1]-geo.begend[0]-1;
      tordata.nunit2 = geo.begend[1]-geo.begend[0];
   }
   tcount = tordata.nunit1+tordata.nunit2; 
   for(i=1;i<tcount;i++) //NEW 18feb99 - set tordata.angsign[i] =0 so will find correct sign
      tordata.angsign[i]=0; 
   count = 0; // counts unit's (1 and 2)
   for(i=1;i<=geo.begend[1]-geo.begend[0];i++){
      if(geo.dir ==0){ // N to c direction
         ii = geo.begend[0]+i-1;
         tordata.cis[i] = geo.cis[ii+1];  // set cis corresponding to Ca
#if 0
if(tordata.cis[i]==1)
  printf(" For i = %d, using cis tordata\n",i);
#endif
//fprintf(fpout,"geo.rtype[%d] = %d\n",ii,geo.rtype[ii]);
         count = count+1;
//         tordata.ttype[count]= geo.rtype[ii];  //THIS IS OLD VERSION (1SEP98) AND (I think) is wrong constaint for bond preceding type 2 unit (NOTE: first unit2 has no constraint - give arbitrary value
         tordata.ttype[count]= geo.rtype[ii+1];  // NEW constaint for bond preceding type 1 (=phi) unit (NOTE: first unit2 has no constraint - give arbitrary value
         count = count+1;
         tordata.ttype[count]= geo.rtype[ii+1];  // NEW constaint for bond preceding type 2 (=psi) unit
         for(k=0;k<=2;k++){
           tordata.unit[1][i][0][k]=geo.pos[ii][2][k]; /* Co pos */
           tordata.unit[1][i][1][k]=geo.pos[ii+1][4][k]; /* N pos of next residue3 */
           tordata.unit[1][i][2][k]=geo.pos[ii][3][k]; /* O pos */
         }             
         if(i<geo.begend[1]-geo.begend[0]){
           for(k=0;k<=2;k++){
             tordata.unit[0][i][0][k]=geo.pos[ii+1][1][k]; /* Ca pos */
             tordata.unit[0][i][1][k]=geo.pos[ii+1][5][k]; /* Cb pos */
           }
         } 
      } // dir = 0
      else {// C to N
         ii =geo.begend[1]-i;
         count = count+1;
         tordata.ttype[count]= geo.rtype[ii];  // constaint for bond preceding type 2 unit
         count = count+1;
         tordata.ttype[count]= geo.rtype[ii];  // constaint for bond preceding type 1unit
         for(k=0;k<=2;k++){
           tordata.unit[1][i][0][k]=geo.pos[ii][2][k]; /* Co pos */
           tordata.unit[1][i][2][k]=geo.pos[ii][3][k]; /* O pos */
           tordata.unit[1][i][1][k]=geo.pos[ii+1][4][k]; /* N pos of next residue3 */
         }             
         if(i<geo.begend[1]-geo.begend[0]){
           for(k=0;k<=2;k++){
             tordata.unit[0][i][0][k]=geo.pos[ii][1][k]; /* Ca pos */
             tordata.unit[0][i][1][k]=geo.pos[ii][5][k]; /* Cb pos */
           }
         } 
      } // dir = 1
   }  // i
   if(geo.dir ==0){// N to c direction
     for(k=0;k<=2;k++){ // add Ca[1] to first and last unit2
        tordata.unit[1][1][3][k]=geo.pos[geo.begend[0]][1][k];
        tordata.unit[1][tordata.nunit2][3][k]=geo.pos[geo.begend[1]][1][k];
     }
   }
   else{ // C to N
      for(k=0;k<=2;k++){ // add Ca[1] to first and last unit2
         tordata.unit[1][1][3][k]=geo.pos[geo.begend[1]][1][k];
         tordata.unit[1][tordata.nunit2][3][k]=geo.pos[geo.begend[0]][1][k];
      }
   }
   if(geo.dir ==0){
      findtorphipsi(tordata);
      findangsign(geo,tordata); // finds sign of torsional force - if needed for phi/psi constraint
   }
   else if(geo.dir ==1){
      findrevtorphipsi(tordata);
      findangsign(geo,tordata); // finds sign of torsional force - if needed for phi/psi constraint
   }
  ns=tordata.sidedata[0][0] = 0; // counts total number of side chain atoms
   for(i=1;i<=tordata.nunit1;i++){ // Run through each unit1, set side chain date
     tordata.sidedata[i][2] = -1; //New required for new copyautotogeo  Indicate geo.iside >0
     tordata.sidedata[i][3]= 1; // NEW 29feb00 - now needed
     if(geo.dir ==0) // N to c direction
         ii = geo.begend[0]+i-1;
     else // C to N
         ii =geo.begend[1]-i-1; // This make ii+1 OK for this direction
//fprintf(fpout,"ii=%d geo.iside[ii+1]=%d\n",ii,geo.iside[ii+1]);
     if(geo.iside[ii+1] ==1){ // if only one side chain atom (= Cb, then no atoms in uside
           tordata.sidedata[i][0] = 0;
     }  
     else{ // geo.iside>1 - there are side chain atoms in uside
        ns=ns+1; 
        tordata.sidedata[i][1]= ns;// start number for this chain
        for(k=0;k<=2;k++)
           tordata.uside[ns][0][k] = geo.pos[ii+1][1][k]; // first atom is Ca
        atomn = 4;  // side chain atoms in geo are numberd 5 (Cb),6,7...
        for(j=1;j<geo.iside[ii+1];j++){ // geo.nside = # of side chain atoms
           atomn = atomn+1;  //  5,6,7 ...
           ns = ns+1; //next atom (Cb, etc.)
           if(j < geo.iside[ii+1] -1){  //check if last two side atoms   (they go on one unit)
              for(k=0;k<=2;k++)
                tordata.uside[ns][0][k] =geo.pos[ii+1][atomn][k];
           } 
           else { // last 2 atoms  
              for(k=0;k<=2;k++){
                tordata.uside[ns][0][k] =geo.pos[ii+1][atomn][k]; 
                tordata.uside[ns][1][k] =geo.pos[ii+1][atomn+1][k]; 
              }
              j = j+1; // Since added two atoms
           }
         }
         tordata.sidedata[i][0] = geo.iside[ii+1]; //e.g. if nside =2; then first side is Ca, second is Cb and Cb2
      } // else
   } //i     
   tordata.sidedata[0][0]= ns; // total # side units
   findsidechi(geo.dir,tordata,geo);
}

int checkanglimit(struct tordyn& tordata)
/* return 1 if phi or psi starts to grow to large, indicates
   that torque consttraint has broken down
 If 1, then quit dynamics
*/
{
   int i,tcount;
   float ang;

   tcount = tordata.nunit1+tordata.nunit2; 
   for(i=1;i<tcount;i++){  
     ang = tordata.phipsi[i];
     if(fabs(ang)>2.0*PI){
//if(verbose>=2)
  fprintf(fpout,"ERROR: bad angle in checkanglimit: res = %d angtype = %d  ang = %5.3f  \n",
        i,tordata.ttype[i],ang);
        return 1;
     }
   }
   return 0;
} 

#if 1 // General Case  GOOD!!! 
void tordyn(struct griddata& den1,struct geometry& geo,
     struct tordyn& tordata)
/*
 New: 17may99 - add case for geo.fixend = 3 - middle of chain
   (see notes at top for details)
   Define new tordata2:
      tordata for Nterm end movement, tordata2 for cterm end movment
 NEW VERSION - if geo.oneside>0   - then just run dynamics for side chain = iside = geo.oneside-1
 begend[0,1] are beginning and end resiudes in geo that
    are run under torsional dynamics (below: ibeg and iend do not = 0)
 geo.begend[0,1] = ends of residues used for tor. dyn.
   ends move depending on valule of fixend
      if bestden does not increase(<1.02*oldbestden), then quit and use previous density
  If fixend = 0, use all atoms, set begend[0,1] = geo.nend[0,1]
  If fixend = 1, fix C term end, input begend[1] = first fixed residue
     and set begend[0] = geo.nend[0]
  If fixend = 2, fix N term end, input begend[0] = last fixed residue
     and set begend[1] = geo.nend[1]
  If fixend = 3  -both ends fixed  - not working yet.
 master routine for Interactive running auto dynamcis.
  NOTE: assumes geo.nend[0] = 1;
Define adata[]  
   adata[1]=  1  first time called, allocate arrays using default NDYN
              2  assign geo to unit, (if # residues < default, if > default- delete and relallocate
              3  just do dyn.
   adata[2]= trace # that are working on
   adata[3]= # steps before returning to main program and redrawing
   adata[4] = # of new units to add to fit.THIS IS NOT USED
   adata[5]= 1 thin extend side, if = 0 do not.
Define  1=ca,2=c,3=o;4=n,5=cb
*/
{
//   static int itry= 0;
   int iside,i,j,k,i1,i2;
   int cenmove;
   float rr;
   static int icall=0;
   int origdiamond;
//   float savpos[NAFIT][NAA][3];

icall = icall+1;
geo.quit = 0; // reset to OK value
if(verbose>=2)
   printf("icall = %d\n",icall);
origdiamond = tordata.usediamond;
//tordata.usediamond = 0;  // if = 1 then use diamond force
//tordata.usediamond = 1;  // if = 1 then use diamond force
if(tordata.usediamond==1){
   if(geo.denmult < -1000.0){ // then must scale density
#if 1 // use if fit is already known
         scaleden(den1,geo,&rr);
#endif
#if 0 // default values determined previously using fit
         scaleden0(den1,geo);
#endif
   }
}
//printf("tordata.adata[1] = %d tordata2.adata[1] = %d\n",tordata.adata[1],tordata2.adata[1]); 
   if(geo.oneside >0){ // just do one side chain dynamics
      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;
   }
   else if (geo.fixend == 1){ // Fix C terminal end
     geo.begend[1] = geo.nend[0]+geo.endmove; // beginning of fixed postiion
     if(geo.begend[1]>=geo.nend[1]){  // If not enough residues, allow all to move
         geo.fixend = 0;
         geo.begend[1]=geo.nend[1];
     }
     geo.dir =1; // C to N direction
     geo.begend[0]=geo.nend[0];
   }
   else if (geo.fixend == 2){ // Fix N terminal end
     geo.begend[0] = geo.nend[1]-geo.endmove-1; // beginning of fixed postiion
     if(geo.begend[0]<=geo.nend[0]){  // If not enough residues, allow all to move
         geo.fixend = 0;
         geo.begend[0]=geo.nend[0];
     }
     geo.dir = 0;// N to C direction
     geo.begend[1]=geo.nend[1];
   }
   else if( geo.fixend ==0){ // let all atoms move
      geo.dir = 0;// N to C direction
      geo.begend[0]=geo.nend[0];
      geo.begend[1] = geo.nend[1];
   }
   else if( geo.fixend ==3) { // Dynamics in center - tordata = nterm, tordata2 = cterm - allocate tordata2 here, tordata below
//     cenmove = (0.6 + geo.endmove/2.0);//  # of res. moved on each side
     cenmove = geo.endmove;//  # of res. moved on each side
#if 0 // This is too stringent 
     if( (geo.nend[0] >= (geo.fnum-1) ) || (geo.nend[1]<= geo.fnum+1) ){
         printf("ERROR: geo.fnum = %d is not in middle\n",geo.fnum);
         goto quitrout;
     }
#endif
     if( (geo.fnum <= geo.nend[0])||( geo.fnum >= geo.nend[1]) ){
         printf("ERROR: geo.fnum = %d is not in middle\n",geo.fnum);
         goto quitrout;
     }
     geo.begend[1] = geo.fnum+cenmove+1; //value Cterm end (tordata2)  beginning of fixed postiion
     if(geo.begend[1]>=geo.nend[1]){  // If not enough residues, allow all to move
         geo.begend[1]=geo.nend[1];
     }
     geo.begend[0]=geo.fnum-1; //reset begend[0] to  value for cterm end
     tordata2.fend = geo.begend[1]-geo.begend[0];
     geo.dir =1; // C to N direction - needed in ncopygeotoauto
     if(tordata2.adata[1]==1){ // allocate tordata2 -use default value for unit,
        if( (geo.begend[1]-geo.begend[0]) > tordata2.maxaunits -2){
          tordata2.maxaunits=geo.begend[1]-geo.begend[0]+2;
          if(tordata2.maxaunits> NUNIT1 -2){
             fprintf(fpout,"ERROR:tordata2.maxaunits = %d > NUNIT1 -2 IN tordyn, Increase NUNIT and recompile\n",
                             tordata2.maxaunits);
             exit(0);
          } 
          tordata2.maxsunits=5*tordata2.maxaunits;
        }
#if 0 // Only tordata is allocated in setupauto -WRONG (18MAY00) BOTH ARE ALLOCATED
        allocatetordata(tordata2);
#endif
        tordata2.adata[1] =2;
     } // tordata2.adata = 1
     if(tordata2.adata[1]==2){
        sidenumreg(geo.begend[0]+1,geo.begend[1],den1,geo); //(add 1 to end[0] because that side chain belongs to ntermend finds sidenum 
        ncopygeotoauto(geo,tordata2); 
        if(geo.quit==1){
           if(verbose>=1)
              printf("geo.quit = 1 after ncopygeotoauto, quiting tordyn\n");
           goto quitrout;
        }
        for(i=0;i<=tordata2.nunit1+1;i++){ // set main chain atom vel = 0
         for(i1=0;i1<=1;i1++)
           for(i2=0;i2<=2;i2++)
             for(k=0;k<=2;k++)
                tordata2.unitv[i1][i][i2][k]=0.0;
        }
        for(i=1;i<=tordata2.nunit1+tordata2.nunit2+1;i++){ /* set initival mainchain Y vector=0 */
          for(k=0;k<=5;k++)
             tordata2.yy[i][k]=0.0; /* set c.m. vel. and ang. rot =0 inititally*/
          tordata2.dq[i]=0.0; /* velocity of rotation about bond from i to j*/   
        }
#if 0 // ERROR: does not set all vel. = 0
        count=0;
        for(i=1;i<=tordata2.nunit1;i++){ 
         for(j=1;j<=tordata2.sidedata[i][0];j++){
            count=count+1;
            for(k=0;k<=5;k++)/* set side chain Y vector = 0 */
              tordata2.yys[count][k]=0.0;
            tordata2.dqs[count]=0;// set dq =0
            for(i1=0;i1<12;i1++)// set side atom vel = 0 
              for(k=0;k<=2;k++)
                tordata2.usidev[count][i1][k]=0.0;
          }
        }
#endif
#if 1 // New (6mar00) version
       for(i=1;i<=tordata2.nunit1;i++){ /* generate rest of dyys and ddqs recursively */
           if(tordata2.sidedata[i][0]>0){
              for(j=tordata2.sidedata[i][1];j<tordata2.sidedata[i][1]+tordata2.sidedata[i][0];j++){/* start at Cb*/
                 for(k=0;k<6;k++)
                    tordata2.yys[j][k]=0.0;
                 tordata2.dqs[j]=0.0;
                 for(i1=0;i1<12;i1++)// set side atom vel = 0 
                   for(k=0;k<=2;k++)
                     tordata2.usidev[j][i1][k]=0.0;
              } // j
              tordata2.omega[i]=0.0;  //set rotation of ring about Cb = 0
          }
       }// i
#endif
        tordata2.adata[1]=3;
      } // tordata2.adata[1]=2
      if(tordata2.adata[1]==3){
        geo.begend[0] = geo.fnum-cenmove-1; // Define geo for nterm half - used to allocate initiliixw tordata below
        if(geo.begend[0]<=geo.nend[0]){  // If not enough residues, allow all to move
            geo.begend[0]=geo.nend[0];
        }
        geo.begend[1]=geo.fnum;
        tordata.fend = geo.begend[1]-geo.begend[0];
        geo.dir=0;
      }
   } // fixend = 3
//goto quitrout;
if(verbose>=2)
  fprintf(fpout,"Runing Interactive Torsional Dyn:icall =%d:fixend = %d  begend = %d %d dir = %d rfor = %5.3f endmove = %d\n",
        icall,geo.fixend,geo.begend[0],geo.begend[1],geo.dir,tordata.rfor,geo.endmove);
   if(tordata.adata[1]==1){ // allocate arrays -use default value for unit,
     if(verbose>=2) 
         fprintf(fpout,"Allocating dengr array dimensions %d %d %d\n",
                           den1.griddim[0],den1.griddim[1],den1.griddim[2]);
      scmat4(den1.dengr,den1.griddim[0],den1.griddim[1],den1.griddim[2],3);// allocates array 3 dim array (mat)
                     // of type  unsigned char - n1,n2,n3,n4 size of each dimension
      finddengrad(den1);
#if 0 // Now done for tordata.adata[1]=2
      if( (geo.begend[1]-geo.begend[0]) > tordata.maxaunits -2){
          tordata.maxaunits=geo.begend[1]-geo.begend[0]+2;
          if(tordata.maxaunits> NUNIT1 -2){
             fprintf(fpout,"ERROR:tordata.maxaunits = %d > NUNIT1 -2 IN tordyn, Increase NUNIT and recompile\n",
                             tordata.maxaunits);
             exit(0);
          } 
          tordata.maxsunits=5*tordata.maxaunits;
      }
#endif
#if 0 // Now done in setupauto
      allocatetordata(tordata);
#endif
      tordata.adata[1] =2;
   }
   if(tordata.adata[1]==2){
#if 1 //initialize tordata    Now only find sidenum on intitial run
      if(geo.fixend!=3){
         sidenumreg(geo.begend[0],geo.begend[1],den1,geo); // finds sidenum 
       }
       else // fixend = 3, this is for nterm end, do not want side chain for last residue because it belongs to cterm end
         sidenumreg(geo.begend[0],geo.begend[1]-1,den1,geo); // finds sidenum
#endif
       ncopygeotoauto(geo,tordata);
#if 0 // checks geometry when call tordyn -ALL and then exits
if(geo.fixend!=3){
  checkstruc(geo,tordata);
  exit(0);
}
#endif
       if(geo.quit==1){
           if(verbose>=1)
              printf("geo.quit = 1 after ncopygeotoauto, quiting tordyn\n");
           goto quitrout;
       }
       if( (geo.begend[1] - geo.begend[0]) > tordata.maxaunits-2){
          deallocatetordata(tordata);
          tordata.maxaunits=geo.begend[1]- geo.begend[0]+2;
          if(tordata.maxaunits> NUNIT1 -2){
             fprintf(fpout,"ERROR:tordata.maxaunits = %d > NUNIT1 -2 IN tordyn, Increase NUNIT and recompile\n",
                             tordata.maxaunits);
             exit(0);
          } 
          tordata.maxsunits=5*tordata.maxaunits; 
          allocatetordata(tordata);          
       }
       for(i=0;i<=tordata.nunit1+1;i++){ // set main chain atom vel = 0
          for(i1=0;i1<=1;i1++)
            for(i2=0;i2<=2;i2++)
              for(k=0;k<=2;k++)
                tordata.unitv[i1][i][i2][k]=0.0;
       }
       for(i=1;i<=tordata.nunit1+tordata.nunit2+1;i++){ /* set initival mainchain Y vector=0 */
          for(k=0;k<=5;k++)
             tordata.yy[i][k]=0.0; /* set c.m. vel. and ang. rot =0 inititally*/
          tordata.dq[i]=0.0; /* velocity of rotation about bond from i to j*/   
       }
#if 0 // Error does not set all vel. = 0
       count=0;
       for(i=1;i<=tordata.nunit1;i++){ 
         for(j=1;j<=tordata.sidedata[i][0];j++){
            count=count+1;
            for(k=0;k<=5;k++)/* set side chain Y vector = 0 */
              tordata.yys[count][k]=0.0;
            tordata.dqs[count]=0;// set dq =0
            for(i1=0;i1<12;i1++)// set side atom vel = 0 
              for(k=0;k<=2;k++)
                tordata.usidev[count][i1][k]=0.0;
         }
       }//i 
#endif
#if 1 // New (6mar00) version
       for(i=1;i<=tordata.nunit1;i++){ /* generate rest of dyys and ddqs recursively */
           if(tordata.sidedata[i][0]>0){
              for(j=tordata.sidedata[i][1];j<tordata.sidedata[i][1]+tordata.sidedata[i][0];j++){/* start at Cb*/
                 for(k=0;k<6;k++)
                    tordata.yys[j][k]=0.0;
                 tordata.dqs[j]=0.0;
                 for(i1=0;i1<12;i1++)// set side atom vel = 0 
                   for(k=0;k<=2;k++)
                     tordata.usidev[j][i1][k]=0.0;
              } // j
              tordata.omega[i]=0.0;  //set rotation of ring about Cb = 0
          }
       }// i
#endif
       tordata.adata[1]=3;      
    } //tordata.adata[1]=2
#if 0//  // debug = run just one step
tordata.adata[3]=1;
#endif
    if(tordata.adata[1]==3){ /* do torsion dyancmaisc*/
      for(i=1;i<=tordata.adata[3];i++){
        if(geo.oneside ==0){
           if(geo.fixend !=3){ // not runing in middle
              runtordyn(den1,geo,tordata,tordata2);
              if(geo.quit==1){
                printf(" geo.quit =1 after runtordyn for geo.fixend!=3\n");
                goto quitrout;
              }
           }
           else{ // geo.fixend =3 alternate runing both ends
             tordata.fixend= tordata2.fixend =3; // This turns on force to keep ends close
             geo.fixend= 2;
             geo.begend[0] = geo.fnum-cenmove-1; // Define geo for nterm half - used to allocate initiliixw tordata below
             if(geo.begend[0]<=geo.nend[0]){  // If not enough residues, allow all to move
                 geo.begend[0]=geo.nend[0];
                 geo.fixend=0;
              }
              geo.begend[1]=geo.fnum;
              geo.dir = 0;// N to C direction
//printf(" Run fuse  Dyn: geo.dir = 0  geo.fnum = %d cenmove = %d\n",geo.fnum,cenmove);
              runtordyn(den1,geo,tordata,tordata2);
              if(geo.quit==1){
                printf(" geo.quit =1 after runtordyn for geo.fixend = 3, tordata\n");
//exit(0);
                 geo.fixend=3; // reset to original value
                 tordata.fixend = tordata2.fixend = 0; // reset to default value
                goto quitrout;
              }
              if(i== tordata.adata[3]){ // then save geo and fpos (position of end)
                 copyautotogeo(geo,tordata);//I think this is redundant if side is a.a. residue
                 copypostofpos(geo,geo.fnum,0);// "0 - for copying the nterm end
              }
              geo.fixend = 1; // geo.fixend for cterm end moves (nterm fixed)
              geo.begend[1] = geo.fnum+cenmove+1; //Cterm end (tordata2)  beginning of fixed postiion
              if(geo.begend[1]>=geo.nend[1]){  // If not enough residues, allow all to move
                 geo.begend[1]=geo.nend[1];
                 geo.fixend=0;
              }
              geo.begend[0]=geo.fnum-1;
              geo.dir =1; // C to N direction
              j=1;
//printf(" Run fuse  Dyn: geo.dir = 1\n");
              runtordyn(den1,geo,tordata2,tordata);
              if(geo.quit==1){
                printf(" geo.quit =1 after runtordyn for geo.fixend = 3, tordata2\n");
//exit(0);
                 geo.fixend=3; // reset to original value
                 tordata.fixend = tordata2.fixend = 0; // reset to default value
                goto quitrout;
              }
              if(i== tordata.adata[3]){ // then save geo and fpos (position of end)
                 copyautotogeo(geo,tordata2);
                 copypostofpos(geo,geo.fnum,1);// "1 - for copying the nterm end
              }
              geo.fixend=3; // reset to original value
              tordata.fixend = tordata2.fixend = 0; // reset to default value
           }// geo.fixend = 3
        } // oneside =0
        else{ // just run 1 side chain dynaimcs
//printf("Running one side dynamics\n");
           runonesidetordyn(1,den1,geo,tordata);
       }
       if(checkanglimit(tordata) ==1){// If bad angle - try again wiht original pos
         if(verbose>=1)
            printf("ERROR  AT 5  i = %d checkanglimit =1  setting geo.quit = 1\n",i);
          geo.quit =1;
//exit(0);
       }
     }// i
     if(geo.quit==1){
           printf("ERROR:  geo.quit = 1 after tordyn - quitting\n");
           goto quitrout;
     }
     if(geo.quit==1){
          fprintf(fpout,"QUITING DYNAMICS BECAUSE  geo.quit = 1\n");
          goto quitrout;
     }
     if(geo.fixend!=3)
          copyautotogeo(geo,tordata);
#if 0// debug  - stop after one round of callling dynamcis
   tordata.adata[1]=4;  // debug  - stop after one step
#endif
#if 0
float phi,psi,lastden;
int checkset = -1; // do not check set
//int nbad;
int rtype = geo.rtype[ geo.nend[1]-1 ];
findphipsi(geo.pos,geo.nend[1]-2,&phi,&psi);
checkresidue(checkset,geo,den1,geo.set,lastden,nbad,phi,psi,rtype); // then end satifies end contitions
#endif 
#if 1  // debug
float avecb2,lastden;
int nbad;
#if 0 //ERROR - SHOULD BE geo.begend???  Used up to 25sep00
   avedenreg(den1,geo,geo.pos,&lastden,&avecb2,&nbad,
                  geo.nend[0],geo.nend[1]); // Find average density for last cterm residue additions
   if(verbose>=2)
     printf(" Length =%d   finalden = %5.3f avecb2 = %5.3f nbad = %d \n",
      geo.nend[1]-geo.nend[0],lastden,avecb2,nbad);
#endif
   avedenreg(den1,geo,geo.pos,&lastden,&avecb2,&nbad,
                  geo.begend[0],geo.begend[1]); // Find average density for last cterm residue additions
   if(verbose>=2)
//   if(verbose>=1)
     printf(" Length =%d   finalden = %5.3f avecb2 = %5.3f nbad = %d \n",
      geo.begend[1]-geo.begend[0],lastden,avecb2,nbad);
#endif 
  } // adata =3
  if(tordata.adata[1]==4) // Used only for debugging, keeps routine fron doeing anythin
    goto quitrout;
  quitrout:;  
//exit(0);
  tordata.usediamond = origdiamond;
}
#endif

void tordynbat(int quitcond,struct griddata& den1,struct geometry& geo,
     struct tordyn& tordata,float *returnden,int *returnbad)
/*
 NOTE:  The way sidenum is handled:
    after each set of steps, call sidenumreg- which is useed to sidechain density,
     But do not run copygeotoauto except for initial call, thus, no not change
   number of sidechain atoms turn batch run
  dir = 0 0r 1:  N to C or C to N direction
  quitcont = 0 (use entire density); 1 (n term);2(c term) as condition for quitting
    quitcond >2 (say) 40 then run 40 rounds of rfinement and use last pos
 master routine for Non-Interactive running auto dynamcis.
   Repeated until reach quit condition:
     1) minimum 5 cycles;
     2) if ((aveden - oldaveden)/fabs(aveden))< .01)
          quits if density change < 1% or density decreases > 1%
   Differs from interactive:
      As in interactive, if adata[1] = 1 then allocate and set adata[1] = 2;
      But, in this case, always assumed new peptide  (that is use step
       that corresponds to adata[1] = 2 in interactive)
  NOTE: assumes geo.nend[0] = 1;
Define adata[]  
   adata[1]=  1  first time called, allocate arrays using default NDYN
              2  assign geo to unit, (if # residues < default, if > default- delete and relallocate
              3  just do dyn.
   adata[2]= trace # that are working on
   adata[3]= # steps before stopping to check density
   adata[4] = # of new units to add to fit. THIS IS NOT USED
   adata[5]= 1 thin extend side, if = 0 do not.
 New condition for stopping:
    If aveden < bestden
     Then run 5 more times,
 geo.begend[0,1] = ends of residues used for tor. dyn.
   ends move depending on valule of fixend
      if bestden does not increase(<1.02*oldbestden), then quit and use previous density
  If fixend = 0, use all atoms, set begend[0,1] = geo.nend[0,1]
  If fixend = 1, fix C term end, input begend[1] = first fixed residue
     and set begend[0] = geo.nend[0]
  If fixend = 2, fix N term end, input begend[0] = last fixed residue
     and set begend[1] = geo.nend[1]
  If fixend = 3  -both ends fixed  - not working yet.
 */
{
   int i,j,k,i1,i2,count,nbad,ncycle,repeatcycle,iside;
   float aveden,bestden,oldbestden,testden,avecb2,bestbad;
   float savpos[NAFIT][NAA][3];
   float rotsav[NAA][3];
   int origfixend;
   int nlim;
static int icall=0; // Used for debugging, keeps track of call number tothis routine
int useold;

useold = 0; // if 1 use old version
geo.quit = 0; // reset to OK value
//tordata.gclin=0.0;
icall = icall+1;
   origfixend = geo.fixend;
float rr;
//tordata.usediamond = 0;  // if = 1 then use diamond force
//tordata.usediamond = 1;  // if = 1 then use diamond force
if(tordata.usediamond==1){
   if(geo.denmult < -1000.0){ // then must scale density
#if 1 // use if fit is already known
         scaleden(den1,geo,&rr);
#endif
#if 0 // default values determined previously using fit
         scaleden0(den1,geo);
#endif
   }
}
   copypos(geo.iside,geo.nend[1]-geo.nend[0],geo.nend[0],geo.nend[0],geo.pos,savpos);// Note: same start point so can keep track of geo.iside[],  initial savpos, if den increases, get new value
   if(geo.oneside >0){ // just do one side chain dynamics
      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;
   }
   else if (geo.fixend == 1){ // Fix C terminal end
     geo.begend[1] = geo.nend[0]+geo.endmove; // beginning of fixed postiion
     if(geo.begend[1]>=geo.nend[1]){  // If not enough residues, allow all to move
         geo.fixend = 0;
         geo.begend[1]=geo.nend[1];
     }
     geo.dir =1; // C to N direction
     geo.begend[0]=geo.nend[0];
   }
   else if (geo.fixend == 2){ // Fix N terminal end
     geo.begend[0] = geo.nend[1]-geo.endmove-1; // beginning of fixed postiion
     if(geo.begend[0]<=geo.nend[0]){  // If not enough residues, allow all to move
         geo.fixend = 0;
         geo.begend[0]=geo.nend[0];
     }
     geo.dir = 0;// N to C direction
     geo.begend[1]=geo.nend[1];
   }
   else if( geo.fixend ==0){ // let all atoms move
      geo.dir = 0;// N to C direction
      geo.begend[0]=geo.nend[0];
      geo.begend[1] = geo.nend[1];
   }
if(verbose>=2)
  fprintf(fpout,"Runing BATCH Torsional Dyn:icall =%d:fixend = %d  begend = %d %d dir = %d quitcond = %d\n",
        icall,geo.fixend,geo.begend[0],geo.begend[1],geo.dir,quitcond);
  bestden = - 1000.0;
//fprintf(fpout,"geo.nend[0]= %d geo.nend[1]= %d\n",geo.nend[0],geo.nend[1]);
   if(tordata.adata[1]==1){ // allocate arrays -use default value for unit,
     if(verbose>=2)
      fprintf(fpout,"Allocating dengr array dimensions %d %d %d\n",
                           den1.griddim[0],den1.griddim[1],den1.griddim[2]);
      scmat4(den1.dengr,den1.griddim[0],den1.griddim[1],den1.griddim[2],3);// allocates array 3 dim array (mat)
                     // of type  unsigned char - n1,n2,n3,n4 size of each dimension
      finddengrad(den1);
      if( (geo.begend[1]-geo.begend[0]) > tordata.maxaunits-2){
          tordata.maxaunits= geo.begend[1]-geo.begend[0]+2;
          tordata.maxsunits=5*tordata.maxaunits;
      }
      if(tordata.maxaunits> NUNIT1 -2){
             fprintf(fpout,"ERROR:tordata.maxaunits = %d > NUNIT1 -2 IN tordyn, Increase NUNIT and recompile\n",
                             tordata.maxaunits);
             exit(0);
      }
#if 0 // now done initially in setup auto 
      allocatetordata(tordata);
#endif
      tordata.adata[1] =2;
   }
   else if(geo.begend[1] - geo.begend[0] > tordata.maxaunits-2){ // reallocate if too large
          deallocatetordata(tordata);
          tordata.maxaunits=geo.begend[1]- geo.begend[0]+2;
          if(tordata.maxaunits> NUNIT1 -2){
             fprintf(fpout,"ERROR:tordata.maxaunits = %d > NUNIT1 -2 IN tordyn, Increase NUNIT and recompile\n",
                             tordata.maxaunits);
             exit(0);
          } 
          tordata.maxsunits=5*tordata.maxaunits; 
          allocatetordata(tordata);          
   }
#if 0
if(useold!=1)
   ncopygeotoauto(geo,tordata);
else
#endif
   ncopygeotoauto(geo,tordata);// New 20mar00 - needed if a.a. side chian
   if(geo.quit==1){
           if(verbose>=1)
              printf("geo.quit = 1 after ncopygeotoauto, quiting tordyn\n");
           goto quitrout;
    }
//   copygeotoauto(geo,tordata);
    for(i=1;i<=tordata.nunit1+tordata.nunit2+1;i++){ /* set initival mainchain vel=0 */
          for(k=0;k<=5;k++)
             tordata.yy[i][k]=0.0; /* set c.m. vel. and ang. rot =0 inititally*/
          tordata.dq[i]=0.0; /* velocity of rotation about bond from i to j*/   
   }
   count=0;
  repeatcycle = 0;
   ncycle = 0;
   oldbestden = -1000.0;
if(useold==1){ // old velocity set
#if 0 // On 6mar00 - found that this did not set all vel = 0
   for(i=1;i<=tordata.nunit1;i++){ /* set side atom vel = 0 */
         for(j=1;j<=tordata.sidedata[i][0];j++){
            count=count+1;
            for(k=0;k<=5;k++)
              tordata.yys[count][k]=0.0;
            tordata.dqs[count]=0;
         }
   }
#endif
#if 1 // New (6mar00) version
   for(i=1;i<=tordata.nunit1;i++){ /* generate rest of dyys and ddqs recursively */
     if(tordata.sidedata[i][0]>0){
        for(j=tordata.sidedata[i][1];j<tordata.sidedata[i][1]+tordata.sidedata[i][0];j++){/* start at Cb*/
            for(k=0;k<6;k++)
                tordata.yys[j][k]=0.0;
            tordata.dqs[j]=0.0;
        }
        tordata.omega[i]=0.0;  //set rotation of ring about Cb = 0
     }
   }
#endif
}
else{ // use new velocity set
   for(i=0;i<=tordata.nunit1+1;i++){ // set main chain atom vel = 0
         for(i1=0;i1<=1;i1++)
           for(i2=0;i2<=3;i2++)
             for(k=0;k<=2;k++)
                tordata.unitv[i1][i][i2][k]=0.0;
   }
#if 0 // On 6mar00 - found that this did not set all vel = 0
   for(i=1;i<=tordata.nunit1;i++){ /* set side atom vel = 0 */
       for(j=1;j<=tordata.sidedata[i][0];j++){
            count=count+1;
            for(k=0;k<=5;k++)
              tordata.yys[count][k]=0.0;
            tordata.dqs[count]=0;
            for(i1=0;i1<12;i1++)// set side atom vel = 0 
              for(k=0;k<=2;k++)
                tordata.usidev[count][i1][k]=0.0;
       }
   }
#endif
#if 1 // New (6mar00) version
   for(i=1;i<=tordata.nunit1;i++){ /* generate rest of dyys and ddqs recursively */
     if(tordata.sidedata[i][0]>0){
        for(j=tordata.sidedata[i][1];j<tordata.sidedata[i][1]+tordata.sidedata[i][0];j++){/* start at Cb*/
            for(k=0;k<6;k++)
                tordata.yys[j][k]=0.0;
            tordata.dqs[j]=0.0;
            for(i1=0;i1<12;i1++)// set side atom vel = 0 
              for(k=0;k<=2;k++)
                tordata.usidev[j][i1][k]=0.0;
        } // j
        tordata.omega[i]=0.0;  //set rotation of ring about Cb = 0
     }
   }// i
#endif
}
#if 0// Old quitcond = 5 now means run 5 times
   if( (quitcond==5)&&(geo.oneside>0) ){// then do not run dynamics
       avedenoneside(den1,geo,tordata,1,&bestden,&nbad); // find density for side chain = 1
       goto quitwhile;
   }
#endif
   if(quitcond >2){ // then just run n = quitcond cycles and quit
     for(j=1;j<=quitcond;j++){
       for(i=1;i<=tordata.adata[3];i++){
//         geo.quit = 0;  // used as marker for point out of range in findforden
         if(geo.oneside ==0){
           runtordyn(den1,geo,tordata,tordata2);

#if 1
           if(geo.quit==1){
               if(verbose>=2)// 30nov00 - do not print this line
                   printf(" For quitcond >2: NOTE NOTE!!  geo.quit = 1 after tordyn (near boundary?) - quitting\n");
               goto quit1;
           }
           if(checkanglimit(tordata) ==1){
             if(verbose>=1)
                 printf("ERROR: AT 2: checkanglimit = 1  SET geo.quit = 1 and quitting\n");
             geo.quit =1;
             goto quit1;
           }
#endif
         }
         else // just run 1 side chain dynaimcs - do not check ang. limit
            runonesidetordyn(1,den1,geo,tordata); // run jsut for sidenum = 1
         if(geo.quit==1){
if(verbose>=2)
  fprintf(fpout,"IN tordynbat - quit = 1, i = %d  bestden = %f\n",i,bestden);
            goto  quit1;
         }
       }//i
#if 0 // for debugging
       copyautotogeo(geo,tordata);
       avedenreg(den1,geo,geo.pos,returnden,&avecb2,returnbad,geo.begend[0],
                  geo.begend[1]); // just applied to main chain atoms
       printf("aveden = %f  nbad = %d\n",*returnden,*returnbad);
#endif
     } // j
     copyautotogeo(geo,tordata);
//     avedenreg(den1,geo,geo.pos,returnden,&avecb2,returnbad,geo.nend[0],
//                  geo.nend[1]); // ERROR??? - SHOULD BE geo.begend, used up to 24sep00
     avedenreg(den1,geo,geo.pos,returnden,&avecb2,returnbad,geo.begend[0],
                  geo.begend[1]); // just applied to main chain atoms
     quit1:; // go here if out of range - thus do not copy to geo and use oringal onre
   }// quitcont >5
   else{ // quitcont <5
#if 1 // New aug24,99 -move so only calculate firsttime  use initial position as initial bestden
      if(quitcond == 0){
           aveden0(den1,geo,geo.pos,&aveden); // Find average density (ave) for this orientation               
      }
      else if(quitcond==1){
          aventermden(1,den1,geo,geo.nend[0],geo.pos,&aveden,&nbad); //1 = do  use Cb ave nterm density 
      }  
      else if(quitcond==2){
         avectermden(1,den1,geo,geo.nend[1]-1,geo.pos,&aveden,&nbad); // ave cterm density 
         if(geo.checkmaindist ==1){ //NOTE -not used anymore - routine maxfitdistgeo probablly not correct now set aveden to weighted den + nbad + nlim
              maxfitdistgeo(geo,&nlim); // nlim = number of Ca and N more than 1.0 Ang. from main
              aveden = aveden*(1.0 - (float)nbad/geo.torend -(float) nlim/3.0);
         }
      }
     bestden = aveden;
     bestbad = nbad;
     copypos(geo.iside,geo.nend[1]-geo.nend[0],geo.nend[0],geo.nend[0],geo.pos,savpos);//Note: same start point so can keep track of geo.iside[]   save best  geo.pos
#endif
    for(j=1;j<=tordata.nprelim;j++){// Run tordata.nprelim cycles before looking for minimum
//      geo.quit = 0;  // used as marker for point out of range in findforden
      for(i=1;i<=tordata.adata[3];i++){ 
         if(geo.oneside ==0)
           runtordyn(den1,geo,tordata,tordata2);
         else // just run 1 side chain dynaimcs - do not check ang. limit
            runonesidetordyn(1,den1,geo,tordata); // run jsut for sidenum = 1
         if(geo.quit== 1){
if(verbose>=2)
  fprintf(fpout,"IN tordynbat at 1 - quit = 1, i = %d  bestden = %f\n",i,bestden);
if( fabs(bestden + 800.9165) < 0.1)
  printf("\n");
            goto quitwhile;
         }
      }
    }
    while(1){ // Use conditional to jump out of loop
      ncycle = ncycle +1;
      for(i=1;i<=tordata.adata[3];i++){
//         geo.quit = 0;  // used as marker for point out of range in findforden
         if(geo.oneside ==0){
           runtordyn(den1,geo,tordata,tordata2);
#if 1
           if(geo.quit==1){
               printf("NOTE  NOTE!!:  geo.quit = 1 after tordyn at 2 (near boundary?) - quitting\n");
               goto quitwhile;
           }
           if(checkanglimit(tordata) ==1){
             if(verbose>=1)
                 printf("ERROR AT 3, checkanglimit=1  setting geo.quit = 1\n");
              geo.quit =1;
           }
#endif
         }
         else // geo.oneside>0 just run 1 side chain dynaimcs - do not check ang. limit
            runonesidetordyn(1,den1,geo,tordata); // run jsut for sidenum = 1
         if(geo.quit==1){
if(verbose>=2)
  fprintf(fpout,"IN tordynbat - quit = 1, i = %d  bestden = %f\n",i,bestden);
            goto quitwhile;
         }
      }
      if(geo.oneside>0){
         avedenoneside(den1,geo,tordata,1,&aveden,&nbad); // find density for side chain = 1
         if(aveden>bestden){
            bestden = aveden;
            bestbad = nbad;
            copyautotogeo(geo,tordata);
//printf("aveden = %5.3f geopos6(x) = geo.pos[geo.nend[0]+oneside][6][0] = %5.3f\n",aveden,geo.pos[geo.nend[0]+geo.oneside][6][0]);
            copygeotorotsav(geo,geo.oneside,rotsav); //geo.oneside = resn
         }
      } //i
      else{ //geo.oneside = 0 standard procedure, not amino acid side
        copyautotogeo(geo,tordata);
        if(quitcond == 0){
           aveden0(den1,geo,geo.pos,&aveden); // Find average density (ave) for this orientation               
        }
        else if(quitcond==1){
           aventermden(1,den1,geo,geo.nend[0],geo.pos,&aveden,&nbad); //1 = do  use Cb ave nterm density 
        }  
        else if(quitcond==2){
           avectermden(1,den1,geo,geo.nend[1]-1,geo.pos,&aveden,&nbad); // ave cterm density 
           if(geo.checkmaindist ==1){ // set aveden to weighted den + nbad + nlim
              maxfitdistgeo(geo,&nlim); // nlim = number of Ca and N more than 1.0 Ang. from main
              aveden = aveden*(1.0 - (float)nbad/geo.torend -(float) nlim/3.0);
           }
        }
//printf("aveden = %5.3f\n",aveden);
        if(aveden>bestden){
           bestden = aveden;
           bestbad = nbad;
           copypos(geo.iside,geo.nend[1]-geo.nend[0],geo.nend[0],geo.nend[0],geo.pos,savpos);//Note: same start point so can keep track of geo.iside[]   save best  geo.pos
        }
      } // geo.oneside = 0
      if( (ncycle>=5)&&( ((aveden - bestden))< 0.0) ){ // condition for stoping
if(verbose>=2)
  fprintf(fpout,"aveden = %f   bestden = %f\n",aveden,bestden);
         if(oldbestden <0)
            testden = oldbestden/1.02;
         else 
            testden = 1.02*oldbestden;
         if( (repeatcycle ==1)&&(bestden<=testden) )
             goto quitwhile;
         else{
           oldbestden = bestden;
           ncycle=0;
           repeatcycle =1;
         }
      }
   } // while
   quitwhile:;
   if(geo.oneside == 0)
      copypos(geo.iside,geo.nend[1]-geo.nend[0],geo.nend[0],geo.nend[0],savpos,geo.pos);//Note: same start point so can keep track of geo.iside[]   save best  geo.pos 
   else{ // side is a residue
       copyrotsavtogeo(geo,geo.oneside,rotsav); //geo.oneside = resn
//printf("oneside = %d geopos6(x) = geo.pos[geo.nend[0]+oneside][6][0] = %5.3f\n",geo.oneside,geo.pos[geo.nend[0]+geo.oneside][6][0]);
   }
   *returnden = bestden;  // return best cterm density
   *returnbad = bestbad;
 }// quitcont <5
 quitrout:;
 geo.fixend = origfixend; // set back to original value
}

