/*
Have placed in this file - routines from other files (used in batchmode)
   that use screen functions:
     1)savepick(struct geometry& geo
 peptide.c   Routines to create and manipulate peptides
NOTE:  HAVE MADE MAJOR CHANGE FROM OLD MAID PROGRAM:
   pos arrays are now of form pos[residue #][atom type]
     (instead of pos[atom type][residue #])
*/
//#include "drawwindow.h"
//#include "screen.h"
#include "math.h"  
#include "peptide.h"
#include "trace.h"
extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version



#define head2 float m[3][3],float det,float aa[2],float bb[2],float u1[3],float u2[3],float v1[3],float v2[3],float dp[3],float cosang,float u1du1,float u1dv1,float v1dv1,float u2du2,float u2dv2,float v2dv2,float u1du2,float u1dv2,float v1du2,float v1dv2
#define func2 m,det,aa,bb,u1,u2,v1,v2,dp,cosang,u1du1,u1dv1,v1dv1,u2du2,u2dv2,v2dv2,u1du2,u1dv2,v1du2,v1dv2
#define head2n float cc[4],int t[3],int ix,float r[5][3],float m[3][3],float det,float dp[3],float cosang,float u1du1,float u1dv1,float v1dv1,float u2du2,float u2dv2,float v2dv2,float u1du2,float u1dv2,float v1du2,float v1dv2
#define func2n cc,t,ix,r,m,det,dp,cosang,u1du1,u1dv1,v1dv1,u2du2,u2dv2,v2dv2,u1du2,u1dv2,v1du2,v1dv2


void setup(struct geometry& geo) /* creates array delta and epsilon*/
{
   int i,j,k;
   for(i=0;i<=2;i++) /*initialize = 0*/
     for(j=0;j<=2;j++) {
        geo.de[i][j]=0;
        for(k=0;k<=2;k++)
          geo.ep[i][j][k]=0;
     }  
   geo.de[0][0]=geo.de[1][1]=geo.de[2][2]=1;
   geo.ep[0][1][2]=geo.ep[2][0][1]=geo.ep[1][2][0]=1;
   geo.ep[1][0][2]=geo.ep[0][2][1]=geo.ep[2][1][0]= -1;
   geo.al[0]= -1; /*al and de are used in interpolation (findatomden and findpdbden)*/
   geo.al[1]=0;
   geo.da[0]=0;
   geo.da[1]=1;
}

#if 1
void assignlevconst(struct geometry& geo)
/*
NEW - 8nov00 - added special angle for Ca-Cb-Cb bond for ring residues:
    Prevously used  geo.acc=DEGTORAD * 109.1;
       Looked at ave from shad and fah:
          TYR: 115,0;
          PHE = HIS = 113.8
    Define new:  geo.tyr = 115.0;  geo.phe = 113.8
constant angle and bond lenght values from 
  M. Levitt (1973)
*/
{
   geo.tyr = DEGTORAD * 116.0; // Ca-Cb-Cg angle for tyrosine
   geo.phe = DEGTORAD * 113.8;
   geo.cnc=DEGTORAD * 112.6; /* angles */
   geo.anc=DEGTORAD * 123.1;
   geo.bnb=DEGTORAD * 109.1;
   geo.bng=DEGTORAD * 107.7;
   geo.ccn=DEGTORAD * 110.2;
   geo.ccc=DEGTORAD * 112.0;
   geo.acn=DEGTORAD * 112.9;
   geo.acc=DEGTORAD * 109.1;
//geo.phe =geo.tyr= geo.acc; // Set back to old values
   geo.bcc=DEGTORAD * 110.1;
   geo.ccs=DEGTORAD * 115.3;
   geo.ccv=DEGTORAD * 109.3;
   geo.oao=DEGTORAD * 122.2;
   geo.nao=DEGTORAD * 124.3;
   geo.nan=DEGTORAD * 120.0;
   geo.cao=DEGTORAD * 120.4;
   geo.can=DEGTORAD * 115.1;
   geo.aaa=DEGTORAD * 121.2;
   geo.aab=DEGTORAD * 120.2;
   geo.aav=DEGTORAD * 118.5;
   geo.aag=DEGTORAD * 114.6;
   geo.nbn=DEGTORAD * 108.7;
   geo.cbn=DEGTORAD * 122.2;
   geo.abc=DEGTORAD * 120.2;
   geo.aba=DEGTORAD * 119.6;
   geo.bbn=DEGTORAD * 108.9;
   geo.bbc=DEGTORAD * 129.6;
   geo.gbc=DEGTORAD * 125.6;
   geo.bbg=DEGTORAD * 105.1;
   geo.csc=DEGTORAD * 99.1;
   geo.css=DEGTORAD * 105.2;
   geo.agn=DEGTORAD * 126.1;
   geo.agb=DEGTORAD * 132.2;
   geo.ggn=DEGTORAD * 108.7;
   geo.agg=DEGTORAD * 122.2;
   geo.bgg=DEGTORAD * 108.7;
   geo.bbb=DEGTORAD * 108; /* 5 membered ring */
   geo.aaa=DEGTORAD * 120; /* 6 membered ring */   
   geo.cn=1.469; /* bondlengths */
   geo.an=geo.en=1.319;
   geo.bn=1.408;
   geo.gn=1.398;
   geo.ao=geo.eo=1.256;
   geo.cv=1.424;
   geo.av=1.379;
   geo.cc=1.527;
   geo.ac=geo.ec=1.527;
   geo.bc=1.543;
   geo.aa=1.394;
   geo.ab=1.375;
   geo.bb=1.343;
   geo.ag=1.387;
   geo.bg=1.398;
   geo.gg=1.425;
   geo.cs=1.806;
   geo.ss=2.086;
   geo.rotcol[1][0]=geo.rotcol[1][1]=geo.rotcol[1][2]=255;
   geo.rotcol[2][0]=255;geo.rotcol[2][1]=geo.rotcol[2][2]=0;
   geo.rotcol[3][0]=geo.rotcol[3][2]=0;geo.rotcol[3][1]=255;
   geo.rotcol[4][0]=geo.rotcol[4][1]=0;geo.rotcol[4][2]=255;
   geo.rotcol[5][0]=geo.rotcol[5][1]=255;geo.rotcol[5][2]=0;
   geo.rotcol[6][0]=geo.rotcol[6][2]=255;geo.rotcol[6][1]=0;
   geo.rotcol[7][1]=geo.rotcol[7][2]=255;geo.rotcol[7][0]=0;
   geo.rotcol[8][0]=geo.rotcol[8][1]=geo.rotcol[8][2]=255/2;
   geo.rotcol[9][0]=geo.rotcol[9][1]=255/2;geo.rotcol[9][2]=0;
   geo.rotcol[10][0]=geo.rotcol[10][2]=255/2;geo.rotcol[10][1]=0;
}
#endif

void firstaa2(struct geometry& geo)
/*
New version: Makes aa. (ca1 to ca2) in arbitrary orientation, ca1 at 0,0,0,
   No cb found or needed for first aa
Returns caw, cw, etc.
*/
{
   int k;
   float cn[3],co[3],cac[3],th,cdn,sn,th1;

   geo.cis[1]=0; // a standard trans pro 
   geo.caw[1][0]=geo.caw[1][1]=geo.caw[1][2]=0;
   cac[0]=1;cac[1]=cac[2]=0; /* unit vector (x) in direction of ca to c*/
   cn[0]= -cos(geo.tt[1]);/*unit vector in directio of c to n */
   cn[1]=sin(geo.tt[1]);
   co[0]= -cos(geo.tt[2]);/*unit vector in direction of c to 0 o*/
   co[1]= -sin(geo.tt[2]);
   co[2]=cn[2]=0;
   for(k=0;k<=2;k++){
      geo.cw[1][k]=geo.dsqrt[1]*cac[k];
      geo.nw[2][k]=geo.cw[1][k]+geo.dsqrt[2]*cn[k];
      geo.ow[1][k]=geo.cw[1][k]+geo.dsqrt[3]*co[k];
   }
   th=PI-geo.tt[1];/*angle between cac and cn*/
   cdn=cos(th);
   sn=1-cdn*cdn;
   th1=th-PI+geo.tt[3]; /*angle between cac and n to ca2*/
   planevec(cac,cn,cdn,sn,geo.dsqrt[4],th,th1,geo.nw[2],geo.caw[2]); /*geo.position ca2*/
}

void  printphipsi(struct geometry geo)
// print out phi psi angle constraints
{
   fprintf(fpout,"Helix:  %5.3f < phi < %5.3f    %5.3f < psii < %5.3f\n",
     (geo.phi[1]+geo.phi0[1])/DEGTORAD, (geo.phi[1]+geo.phi1[1])/DEGTORAD,
     (geo.psi[1]+geo.psi0[1])/DEGTORAD, (geo.psi[1]+geo.psi1[1])/DEGTORAD);
   fprintf(fpout,"Sheet:  %5.3f < phi < %5.3f    %5.3f < psii < %5.3f\n",
     (geo.phi[2]+geo.phi0[2])/DEGTORAD, (geo.phi[2]+geo.phi1[2])/DEGTORAD,
     (geo.psi[2]+geo.psi0[2])/DEGTORAD, (geo.psi[2]+geo.psi1[2])/DEGTORAD);
}

void findphipsi(float pos[NAFIT][NAA][3],int i,float *phi,float *psi)
/* 
Modified from same routine in old maid
New version: uses findchi routine
   Uses N, (old version used Cb) Now agrees with Frodo result
  Define pos[i] 1=ca,2=c,3=o;4=n,5=cb; 6=cb2
for i: phi is angle formed by Ci-Ni+1-Cai+1-Ci+1
       psi is angle formed by Ni+1-Cai+1-Ci+1-Ni+2
*/
{
   findchi(pos[i][2],pos[i+1][4],pos[i+1][1],pos[i+1][2],phi);
   findchi(pos[i+1][4],pos[i+1][1],pos[i+1][2],pos[i+2][4],psi);
}

void findphipsifit(struct pepfit& fit,int i,float *phi,float *psi)
/* 
Modified to use fit rather then geo as above
*/
{
   findchi(fit.pos[i][2],fit.pos[i+1][4],fit.pos[i+1][1],fit.pos[i+1][2],phi);
   findchi(fit.pos[i+1][4],fit.pos[i+1][1],fit.pos[i+1][2],fit.pos[i+2][4],psi);
}


void findconstants(struct geometry& geo)
/* finds standard distance and angles for amino acid.  dd[1-11] are the
   squared constraints used in molecular dynamics routine
e=carbonyl c;a=trivalent carbon in 6 membered ring  
r1=ec=c-c; r2=en;r3=eo;r4=cn;
for angles, a and e are equivalent
t1=can; t2=cao;t3=anc;t4=nca;
New value using Procheck values:
Engh & Huber mainchain bond angles
The table below shows the standard main-chain bond angles derived from small-molecule data by Engh & Huber (1991), with standard deviations in the rightmost
column. The atom-labelling follows that used in the X-PLOR dictionary, with additional types, defined by Engh & Huber, asterisked. 
Mainchain bond angles
----------------------------------------------------------------
Angle           | X-PLOR labelling               | Value | sigma 
----------------------------------------------------------------
C-N-Calpha      | C-NH1-CH1E    (except Gly,Pro) | 121.7 | 1.8
                | C-NH1-CH2G*   (Gly)            | 120.6 | 1.7
                | C-N-CH1E      (Pro)            | 122.6 | 5.0
                |                                |       |    
Calpha-C-N      | CH1E-C-NH1    (except Gly,Pro) | 116.2 | 2.0
                | CH2G*-C-NH1   (Gly)            | 116.4 | 2.1
                | CH1E-C-N      (Pro)            | 116.9 | 1.5
                |                                |       |    
Calpha-C-O      | CH1E-C-O      (except Gly)     | 120.8 | 1.7
                | CH2G*-C-O     (Gly)            | 120.8 | 2.1
                |                                |       |    
Cbeta-Calpha-C  | CH3E-CH1E-C   (Ala)            | 110.5 | 1.5
                | CH1E-CH1E-C   (Ile,Thr,Val)    | 109.1 | 2.2
                | CH2E-CH1E-C   (the rest)       | 110.1 | 1.9
                |                                |       |    
N-Calpha-C      | NH1-CH1E-C    (except Gly,Pro) | 111.2 | 2.8
                | NH1-CH2G*-C   (Gly)            | 112.5 | 2.9
                | N-CH1E-C      (Pro)            | 111.8 | 2.5
                |                                |       |    
N-Calpha-Cbeta  | NH1-CH1E-CH3E (Ala)            | 110.4 | 1.5
                | NH1-CH1E-CH1E (Ile,Thr,Val)    | 111.5 | 1.7
                | N-CH1E-CH2E   (Pro)            | 103.0 | 1.1
                | NH1-CH1E-CH2E (the rest)       | 110.5 | 1.7
                |                                |       |    
O-C-N           | O-C-NH1       (except Pro)     | 123.0 | 1.6
                | O-C-N         (Pro)            | 122.0 | 1.4
--------------------------------------------------------------  

Previously used 
traherdy angle. M. Levitt has for t4 112
     where ideal is 109.47:
e=carbonyl c;a=trivalent carbon in 6 membered ring  
r1=ec=c-c; r2=en;r3=eo;r4=cn;
for angles, a and e are equivalent
t1=can; t2=cao;t3=anc;t4=nca;
Type 1 = helix, 2 = sheet, 3  = left helix 4= epsilon 5 = delta
*/
{
   int i,k;
   float r1,r2,r3,r4,t1,t2,t3,t4,t5,uv[3],cv[3],cdist;

#if 0 // angles from M. Levitt
   r1=1.527;r2=1.319;r3=1.256;r4=1.469;
   t1=115.1*DEGTORAD;t2=120.4*DEGTORAD;t3=123.1*DEGTORAD;t4=acos(-1.0/3.0);
   t5=112.0*DEGTORAD; /* C-C-C bond angle  Only used in automaid.c */
#endif
#if 1 // Engh and Huber
   r1=1.525;r2=1.329;r3=1.231;r4=1.458;
   t1=116.2*DEGTORAD;t2=120.8*DEGTORAD;
   t3=121.7*DEGTORAD;
   t4=acos(-1.0/3.0); // =109.47 Use ideal tetrahedral (Huber about 109 - 111)
   t5=112.0*DEGTORAD; /* C-C-C bond angle  Only used in automaid.c to build generic side chain*/
#endif
   geo.dda[1]=r1*r1;
   geo.dda[2]=r2*r2;
   geo.dda[3]=r3*r3;
   geo.dda[4]=r4*r4;
   for(i=1;i<=4;i++)
      geo.dsqrt[i]=sqrt(geo.dda[i]);  /*Needed for firstaa() dsqrt is the standard distance*/
   geo.tt[1]=t1;
   geo.tt[2]=t2;
   geo.tt[3]=t3;
   geo.tt[4]=t4;
   geo.tt[5]=t5;
   firstaa2(geo); /*make aa (caw[1]-cw[1]-ow[1]-nw[2]-caw[2] to caw[2], no cb) with ca1 at 0,0,0 */
   for(k=0;k<=2;k++){
     uv[k]=geo.caw[2][k]-geo.caw[1][k];
     cv[k]=geo.cw[1][k]-geo.caw[1][k];
   }
   geo.cadist=sqrt(dot(uv,uv)); /*this used to be called cadist*/
   geo.cardist = 0.0; // This is the runing distance along atoms from Ca to Ca
   for(k=0;k<=2;k++)
      uv[k]=geo.cw[1][k]-geo.caw[1][k];
   geo.cardist = geo.cardist + sqrt( dot(uv,uv) );
   for(k=0;k<=2;k++)
      uv[k]=geo.nw[2][k]-geo.cw[1][k];
   geo.cardist = geo.cardist + sqrt( dot(uv,uv) );
   for(k=0;k<=2;k++)
      uv[k]=geo.caw[2][k]-geo.nw[2][k];
   geo.cardist = geo.cardist + sqrt( dot(uv,uv) );
   printf("geo.cadist = %5.3f   geo.cardist = %5.3f\n",geo.cadist,geo.cardist);
   cdist=sqrt(dot(cv,cv));
   for(k=0;k<=2;k++){
     uv[k]=uv[k]/geo.cadist;
     cv[k]=cv[k]/cdist;
   }
   geo.coneang1=acos(dot(uv,cv));
   geo.coneang2=PI-geo.coneang1-geo.tt[1]; /*angle between C-N and ca-ca*/
   geo.coneang4=PI-geo.coneang2-geo.tt[3];
    geo.dda[5]=r1*r1+r2*r2-2*r1*r2*cos(t1); /*constrains angles t1*/
    geo.dda[6]= r2*r2+r4*r4-2*r2*r4*cos(t3);
    geo.dda[7]=r1*r1+r3*r3-2*r1*r3*cos(t2);
    geo.dda[8]=r4*r4+r1*r1-2*r4*r1*cos(t4);
    geo.dda[9]=2*r1*r1*(1.0-cos(t4)); /*r8 and r9 fix cb*/
    for(k=0;k<=2;k++) /*find O[i]- Ca[i+1] dist*/
       uv[k]=geo.caw[2][k]-geo.ow[1][k];
    geo.dda[10]=dot(uv,uv); /*O - Ca dist - constrains peptide ring*/
    geo.dda[11]=r2*r2+r3*r3-2*r2*r3*cos(2*PI-t1-t2); /* O[i]-N[i+1] dist, const. ring*/
    for(i=5;i<=11;i++)
      geo.dsqrt[i]=sqrt(geo.dda[i]);  /*dsqrt is the standard distance*/
    geo.hbonddist=2.93; /* Distance between O and N in helix hydcrogen bond*/
/* Following "sheet" parameters are angle constraints used in 
      tordyn, classifyphipsi 
   Different values are assigned in extendfit.c++
New version:  keep track of sign of force:
 NOTE! IMPORTANT! IMPORTANT! - max and min angles must be about 5 degreees from
    + or - PI.
    Procedure assign sign of force based on first excursion from min < max
      region.  This first excursion must be  > -PI and < PI.
  1) ang is absolute - ie can get > 180 or less than -180, 
  2) THIS ASSUMES DOES NOT GO COMPLETELTY AROUND WHICH WOULD SCREW THINGS UP -]
       SINCE TORTIONAL FORCE ADDED - SHOULD PREVENT THIS.
  3) Check is sign from previous step = 0, if not, then keep this sign
       if = 0, then test:  if ang < sheet[][0], then sign = 1
          if ang> sheet[][1], then sign = -1;
  4) mag grows from 0 to 1 as ang deviates from min/max 0 to sheet[][4]
                 *mag = fabs(ang - minang)/delang;
        *mag = fabs(ang - maxang)/delang;

*/
//Type 1 = helix, 2 = sheet, 3  = left helix 4= epsilon 5 = delta
//    6= gammar 7 = gammal 

float delang = 30.0*DEGTORAD; // generic allowed +/- range
#if 1 // Put some limits on the 0 constraint case -  -PI <ang < PI
   geo.sheet[0][0]= -DEGTORAD*180.0; /*phi min*/
   geo.sheet[0][1]= DEGTORAD*180.0;  /* phi max*/
   geo.sheet[0][2]= -DEGTORAD*180.0; /*psi min*/
   geo.sheet[0][3]= DEGTORAD*180.0; /*psi max*/
   geo.sheet[0][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
#endif  
/*type 1: helix parameters*/
   geo.phi[1] = -59.0*DEGTORAD; //  Values used in maid: Ideal helix parameters
   geo.psi[1] = -45.0*DEGTORAD;
#if 1 // Used up to 30may00
   geo.sheet[1][0]= -DEGTORAD*120.0; /*phi min*/
   geo.sheet[1][1]= -DEGTORAD*40.0;  /* phi max*/
   geo.sheet[1][2]= -DEGTORAD*65.0; /*psi min*/
   geo.sheet[1][3]= DEGTORAD*10.0; /*psi max*/
   geo.sheet[1][7] = DEGTORAD*60.0; // Used to 30may00  TRY  Increase range   +/- range 
#endif
#if 0 //NOTE - THIS IS NOW TEMPORARILY SET IN TRACE.C++ Used after  30may00 try smaller allowed rang
   geo.sheet[1][0]= -DEGTORAD*110.0; /*phi min*/
   geo.sheet[1][1]= -DEGTORAD*50.0;  /* phi max*/
   geo.sheet[1][2]= -DEGTORAD*45.0; /*psi min*/
   geo.sheet[1][3]= DEGTORAD*10.0; /*psi max*/
   geo.sheet[1][7] = DEGTORAD*30.0; //New 30may00 - try smaller range +/- range 
#endif
/* Try including correlation of phi/psi:  psim = slope*phi + intercept 
       psim-range < psi < psim+range
*/
   geo.sheet[1][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
   geo.sheet[1][5] = -315.0/210.0; //slope
   geo.sheet[1][6]= 135.0*PI/180.0 -geo.sheet[1][5]*(-PI);//22jun00,- try differnt intercept 
//   geo.sheet[1][5] = -280.0/225.0; //slope  - 26apr99  TRY decrease slope
//   geo.sheet[1][6] = geo.psi[1] - geo.sheet[1][5]*geo.phi[1]; // intercept
/*type 2: sheet parameters*/
    geo.phi[2] = -139.0*DEGTORAD;// Ideal value for antiparrel sheet
    geo.psi[2] = 135.0*DEGTORAD;
#if 1 // Used up to 30may00
   geo.sheet[2][0]= -DEGTORAD*165.0; /*phi min*/
   geo.sheet[2][1]= -DEGTORAD*45.0;  /* phi max*/
//(NOTE! used before 18feb99   geo.sheet[2][2]= DEGTORAD*100.0; /*psi min*/
   geo.sheet[2][2]= DEGTORAD*90.0; /*psi min*/
// (NOTE! used before 18feb99, bad because = PI,  means signphipsi in error)   geo.sheet[2][3]=DEGTORAD*180.0; /*psi max*/
   geo.sheet[2][3]=DEGTORAD*175.0; //psi max - New CORRECT - after 18feb99
   geo.sheet[2][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
#endif
#if 0 // Tried after 30may00 using more restrived angles - POORER RESULT
   geo.sheet[2][0]= -DEGTORAD*150.0; /*phi min*/
   geo.sheet[2][1]= -DEGTORAD*60.0;  /* phi max*/
   geo.sheet[2][2]= DEGTORAD*110.0; /*psi min*/
   geo.sheet[2][3]=DEGTORAD*160.0; //psi max - New CORRECT - after 18feb99
   geo.sheet[2][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
#endif
#if 0 // Tried after 30may00 using less restrived angles - NO BETTER
   geo.sheet[2][0]= -DEGTORAD*170.0; /*phi min*/
   geo.sheet[2][1]= -DEGTORAD*35.0;  /* phi max*/
   geo.sheet[2][2]= DEGTORAD*80.0; /*psi min*/
   geo.sheet[2][3]=DEGTORAD*175.0; //psi max - New CORRECT - after 18feb99
   geo.sheet[2][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
#endif
/*type 3: left helix region parameters*/
    geo.phi[3] = 61.0*DEGTORAD; //   Ideal left helix parameters
    geo.psi[3] = 45.0*DEGTORAD;
   geo.sheet[3][0]= DEGTORAD*45.0; /*phi min*/
   geo.sheet[3][1]= DEGTORAD*90.0;  /* phi max*/
   geo.sheet[3][2]= DEGTORAD*0.0; /*psi min*/
   geo.sheet[3][3]=DEGTORAD*65.0; /*psi max*/
   geo.sheet[3][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
#if 0 // Used before 7oct98
/*type 4: epsilon   - this is made a large - catch all region  used for glycine*/
    geo.phi[4] = 90.0*DEGTORAD; //   Ideal epsilon parameters
    geo.psi[4] = -129.0*DEGTORAD;
   geo.sheet[4][0]= DEGTORAD*20.0; /*phi min*/
   geo.sheet[4][1]= DEGTORAD*165.0;  /* phi max*/
   geo.sheet[4][2]= -DEGTORAD*175.0; /*psi min*/
   geo.sheet[4][3]= -DEGTORAD*50.0; /*psi max*/
   geo.sheet[4][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
#endif
#if 1
/*type 4: epsilon   -this is now limited region - not for glycines*/
//printf("Warning : USING LIMITED EPSILON REGION FOR PHI/PSI\n");
    geo.phi[4] = 55.0*DEGTORAD; //   Ideal epsilon parameters
    geo.psi[4] = -135.0*DEGTORAD;
   geo.sheet[4][0]= DEGTORAD*39.0; /*phi min*/
   geo.sheet[4][1]= DEGTORAD*91.0;  /* phi max*/
   geo.sheet[4][2]= -DEGTORAD*179.5; /*psi min*/
   geo.sheet[4][3]= -DEGTORAD*130.0; /*psi max*/
   geo.sheet[4][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
#endif
//type 5: delta  parameters
    geo.phi[5] = -101.0*DEGTORAD; //   Ideal delta parameters
    geo.psi[5] = 61.0*DEGTORAD;
   geo.sheet[5][0]= geo.phi[5]-delang; /*phi min*/
   geo.sheet[5][1]= geo.phi[5]+delang;  /* phi max*/
   geo.sheet[5][2]= geo.psi[5]-delang; /*psi min*/
   geo.sheet[5][3]= geo.psi[5]+delang; /*psi max*/
   geo.sheet[5][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
//type 6: gammar parameters
    geo.phi[6] = -101.0*DEGTORAD; //   Ideal gammar parameters
    geo.psi[6] = 11.0*DEGTORAD;
   geo.sheet[6][0]= geo.phi[6]-delang; /*phi min*/
   geo.sheet[6][1]= geo.phi[6]+delang;  /* phi max*/
   geo.sheet[6][2]= geo.psi[6]-delang; /*psi min*/
   geo.sheet[6][3]= geo.psi[6]+delang; /*psi max*/
   geo.sheet[6][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
//type 7: gammal parameters
    geo.phi[7] = 91.0*DEGTORAD; //   Ideal gammar parameters
    geo.psi[7] = 1.0*DEGTORAD;
   geo.sheet[7][0]= geo.phi[7]-delang; /*phi min*/
   geo.sheet[7][1]= geo.phi[7]+delang;  /* phi max*/
   geo.sheet[7][2]= geo.psi[7]-delang; /*psi min*/
   geo.sheet[7][3]= geo.psi[7]+delang; /*psi max*/
   geo.sheet[7][4]= DEGTORAD*20.0; //  range of which mag varies from 0 to 1
// type 8: Proline helix region:
   geo.phi[8] = -58.0*DEGTORAD;//   Ideal proline parameters
   geo.psi[8] = -35.0*DEGTORAD;
   geo.sheet[8][0]= -82.0*DEGTORAD; /*phi min*/
   geo.sheet[8][1]= -50.0*DEGTORAD;  /* phi max*/
   geo.sheet[8][2]=  -50.0*DEGTORAD; /*psi min*/
   geo.sheet[8][3]=  -16.0*DEGTORAD; /*psi max*/
   geo.sheet[8][4]= DEGTORAD*10.0; //  range of which mag varies from 0 to 1
// type 9: Proline sheet region:
   geo.phi[9] = -58.0*DEGTORAD;//   Ideal proline parameters
   geo.psi[9] = 144.0*DEGTORAD;
   geo.sheet[9][0]= -82.0*DEGTORAD; /*phi min*/
   geo.sheet[9][1]= -50.0*DEGTORAD;  /* phi max*/
   geo.sheet[9][2]=  115.0*DEGTORAD; /*psi min*/
   geo.sheet[9][3]=  170.0*DEGTORAD; /*psi max*/
   geo.sheet[9][4]= DEGTORAD*10.0; //  range of which mag varies from 0 to 1
//type 10: NEW added 13oct00  Lower right corner of helix region  
   geo.phi[10] = 0.0*DEGTORAD; 
   geo.psi[10] = -60.0*DEGTORAD;
   geo.sheet[10][0]= -DEGTORAD*50.0; /*phi min*/
   geo.sheet[10][1]= DEGTORAD*20.0;  /* phi max*/
   geo.sheet[10][2]= -DEGTORAD*90.0; /*psi min*/
   geo.sheet[10][3]= -DEGTORAD*40.0; /*psi max*/
   geo.sheet[10][4]= DEGTORAD*10.0; //  range of which mag varies from 0 to 1
/* These parameters are used for Powell minimization and for checkangrange
   NOTE: IMPORTANT: ideal angle + or - range must be  between  -180 and + 180
   Use ideal value of angle (= phi[i] or psi[i]) give above
   Range is from:  (e.g.)
        phi[i]+phi0[i]  <  phi[i] + phi1[i]
   Range is (in most cases) larger than above (geo.sheet paramaters) because Powell
     minimization imposes an absolute limit on angle, does not allow it
     to go outside range, wherease, tordyn allows any value of angle,
     but imposes maximum force when goes outside range. 
*/
    geo.phi0[1] = -60*DEGTORAD; //Lowerbounds = phi + phi0 = 110   hellix
    geo.psi0[1] = -40.0*DEGTORAD;
    geo.phi1[1] = 30.0*DEGTORAD; //phi + phi1 = 37  helix
    geo.psi1[1] = (geo.psi[5]-geo.psi[1]-1.0*DEGTORAD)/2.0;// prevents overlap with delta (= 5)
    geo.phi0[2] = -40.0*DEGTORAD;// Lowerbounds - sheet
    geo.psi0[2] = -(geo.psi[2]-geo.psi[5] -1.0*DEGTORAD)/2.0;// prevents overlap with delta (= 5
    geo.phi1[2] = 80.0*DEGTORAD;// Upper bounds - sheet
    geo.psi1[2] = 45.0*DEGTORAD;
    for(i=3;i<=7;i++){
       geo.phi0[i] = -50*DEGTORAD; //Lowerbounds = phi + phi0 = 
       geo.psi0[i] = -50.0*DEGTORAD;
       geo.phi1[i] = -50.0*DEGTORAD; //Upperbound = phi + phi1
       geo.psi1[i] = 50.0*DEGTORAD;
   }
   // Use very large range for epsilon
   geo.phi0[4] = -80*DEGTORAD; //Lowerbounds = phi + phi0 = 
   geo.psi0[4] = -50.0*DEGTORAD;
   geo.phi1[4] = 80.0*DEGTORAD; //Upperbound = phi + phi1
   geo.psi1[4] = 110.0*DEGTORAD;
         
    printphipsi(geo);
    geo.plerr = cos(75*DEGTORAD);  //(plane error (rad))= allowed deviation of cos(ang) from 0 where ang in angle
               //between third vector and vector perp. to first 2 vectors.
// these paramters are used in phiallow() and psiallow: Non glycine residues
   geo.phiallow[0][0] = -140.0*DEGTORAD; // min allowed helix phi
   geo.phiallow[0][1] = -40.0*DEGTORAD; // max allowed  helix phi
   geo.psiallow[0][0] = -70.0*DEGTORAD; // min allowed helix psi
   geo.psiallow[0][1] = 40.0*DEGTORAD; // max allowed helix  phi
   geo.phiallow[1][0] = -170.0*DEGTORAD; // min allowed sheet phi
   geo.phiallow[1][1] = -40.0*DEGTORAD; // max allowed  sheet phi
   geo.psiallow[1][0] = 100.0*DEGTORAD; // min allowed sheet psi
   geo.psiallow[1][1] = 180.0*DEGTORAD; // max allowed sheet  phi
   geo.phiallow[2][0] = geo.phiallow[0][0]; // min allowed delta phi
   geo.phiallow[2][1] = geo.phiallow[0][1]; // max allowed  delta phi
   geo.psiallow[2][0] = geo.psiallow[0][1]; // min allowed delta psi
   geo.psiallow[2][1] = geo.psiallow[1][0]; // max allowed delta  phi
   geo.phiallow[3][0] = 40.0*DEGTORAD; // min allowed lefthelix phi
   geo.phiallow[3][1] = 80.0*DEGTORAD; // max allowed  lefthelix phi
   geo.psiallow[3][0] = -10.0*DEGTORAD; // min allowed lefthelix psi
   geo.psiallow[3][1] = 90.0*DEGTORAD; // max allowed lefthelix  phi
}

#if 0  // This is  old version from maid - Gives idenditcal result to new
//        but (confusingly) use +phi, +psi, but negative vectors
void makerevalanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi) 
/* reverse of makealanine2: 
  Old version from maid
    1) Use c to Ca vector and Ca as origin for psi rotatin
    2) Uses Ca to N vector and Ca as orign for phi rotation
       This choice of vectors is equivalnet to -phi and -psi 
        rotation  - which is what is wanted
   NOTE:  I think that any point on the axis of rotation can
     be used as origin
Given  Ca1,C1 abd N2 find a.a. 0 (N1,C0,O0,Ca0, etc) if aanum=1
Does not assume ideal distances and angles in previous a.a.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2
if ifind=1, do not find o or cb
*/
{
   int i,k;
   float mag1,mag2;
   float ct[3];/*intermed co-ord*/
   float cca[3],cc[3],cn[3],nc[3],cac[3],n[3],c[3],cb[3];
   float rom[3][3],cdn,sn,th,th1;
   
   /*here we set phi=psi=0 and find atoms positons*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      cca[k]=pos[aanum][1][k]-pos[aanum][2][k]; /*the c to ca vector*/
      nc[k]=pos[aanum][2][k]-pos[aanum+1][4][k]; /* the n to c vector*/
   }
   mag1=sqrt(dot(cca,cca));
   mag2=sqrt(dot(nc,nc));
   for(k=0;k<=2;k++){
      cca[k]=cca[k]/mag1;
      nc[k]=nc[k]/mag2;
   }
   rotaxis(geo.de,geo.ep,psi,cca,rom); /*rom is the rotation matrix for rotation of phi about 
       unit vector nca, vector nca must pass through origin*/
   /*find n and c for phi=psi=0;*/
   cdn=dot(cca,nc);
   th=acos(cdn);/*theta= angle between 2 basis vectors n and c >0 (anticlockwise)*/
   sn=1-cdn*cdn; /*= sqr(sin(theta)*/
   th1=th+PI-geo.tt[4];
   planevec(nc,cca,cdn,sn,geo.dsqrt[4],th,th1,pos[aanum][1],n);
                        /*return vector n1 in same plane as n2,c and ca*/  
   th1=th-geo.tt[4]-geo.tt[3];                    
   planevec(nc,cca,cdn,sn,geo.dsqrt[2],th,th1,n,c);
   if(ifind!=1){
     /*for cb need define a third vector to define cbeta which is out of plane*/
     for(k=0;k<=2;k++){
        cn[k]=(n[k]-pos[aanum][1][k])/geo.dsqrt[4]; /* the ca-n vector*/
        cac[k]= -cca[k]; /* the ca-c vector*/
     }
     crossprod(cn,cac,cc); /*cc is the component perpendicula to nca,cac*/
     for(i=0;i<=2;i++){ /* the position of cb*/
       cc[i]=cc[i]/sin(geo.tt[4]);
       cb[i]=geo.dsqrt[1]*(-0.5*cn[i]-0.5*cac[i]+sqrt(2.0/3.0)*cc[i])+pos[aanum][1][i];
     }
   /*rotate cb,n,c by psi about c-ca*/
     tranm(cb,pos[aanum][1]); /* puts cw at origin*/
     multvec(rom,cb,pos[aanum][5]);
     tran(pos[aanum][5],pos[aanum][1]);
   }
   tranm(n,pos[aanum][1]);tranm(c,pos[aanum][1]);
   multvec(rom,n,pos[aanum][4]);tran(pos[aanum][4],pos[aanum][1]);
   multvec(rom,c,ct);tran(ct,pos[aanum][1]);
   /*rotate c, by phi about cn*/
   for(i=0;i<=2;i++)
     cn[i]=(pos[aanum][4][i]-pos[aanum][1][i])/geo.dsqrt[4]; /* this is different than cn above*/
   rotaxis(geo.de,geo.ep,phi,cn,rom);
   tranm(ct,pos[aanum][1]);
   multvec(rom,ct,pos[aanum-1][2]); 
   tran(pos[aanum-1][2],pos[aanum][1]); 
   /* find o and ca1 in same plane as c,ca2, and n*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      cca[k]=pos[aanum][4][k]-pos[aanum][1][k]; /*the ca to n vector*/
      nc[k]=pos[aanum-1][2][k]-pos[aanum][4][k]; /* the n to c vector*/
   }
   mag1=sqrt(dot(cca,cca));
   mag2=sqrt(dot(nc,nc));
   for(k=0;k<=2;k++){
      cca[k]=cca[k]/mag1;
      nc[k]=nc[k]/mag2;
   }
   cdn=dot(cca,nc);
   th=acos(cdn);/*theta= angle between 2 basis vectors n and c >0 (anticlockwise)*/
   cdn=cos(th); /* = cn dot nca = cos(theta)*/
   sn=1-cdn*cdn; /*= sqr(sin(theta)*/
   th1=th-PI+geo.tt[1];
   planevec(cca,nc,cdn,sn,geo.dsqrt[1],th,th1,pos[aanum-1][2],pos[aanum-1][1]);
   if(ifind!=1){
     th1=th+geo.tt[2]+geo.tt[1]-PI;
     planevec(cca,nc,cdn,sn,geo.dsqrt[3],th,th1,pos[aanum-1][2],pos[aanum-1][3]); 
   }     
}
#endif

#if 1  // This is new  - Gives identical results as old 
void makerevalanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi) 
/* reverse of makealanine2: 
 This is new  version - Identical results but explicitly use -psi and -phi for rotatino
   Changes:  1) Use ca to c vector  and Ca as origin for psi rotation
             2) Use n to Ca vector and Ca as origin of phi rotatin
   Use -phi and -psi because phi indicate how to form residue 2 give 1
      Here want to go from residue 2 to residue 1
   NOTE:  Any point on the axis of rotation can
     be used as origin - Thus, could use either Ca or N for
     phi rotation orign
Given  Ca1,C1 abd N2 find a.a. 0 (N1,C0,O0,Ca0, etc) if aanum=1
Does not assume ideal distances and angles in previous a.a.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2
if ifind=1, do not find o or cb
*/
{
   int i,k;
   float mag1,mag2;
   float ct[3];/*intermed co-ord*/
   float cca[3],cc[3],cn[3],nc[3],cac[3],n[3],c[3],cb[3],nca[3];
   float rom[3][3],cdn,sn,th,th1;
   
   geo.cis[aanum-1] = 0;  // Indicate is trans proline
   /*here we set phi=psi=0 and find atoms positons*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      cca[k]=pos[aanum][1][k]-pos[aanum][2][k]; /*the c to ca vector*/
      nc[k]=pos[aanum][2][k]-pos[aanum+1][4][k]; /* the n to c vector*/
   }
   mag1=sqrt(dot(cca,cca));
   mag2=sqrt(dot(nc,nc));
   for(k=0;k<=2;k++){
      cca[k]=cca[k]/mag1;
      cac[k]= -cca[k]; /* the ca-c vector  = psi rotation vector*/
      nc[k]=nc[k]/mag2;
   }
   /*find n and c for phi=psi=0;*/
   cdn=dot(cca,nc);
   th=acos(cdn);/*theta= angle between 2 basis vectors n and c >0 (anticlockwise)*/
   sn=1-cdn*cdn; /*= sqr(sin(theta)*/
   th1=th+PI-geo.tt[4];
   planevec(nc,cca,cdn,sn,geo.dsqrt[4],th,th1,pos[aanum][1],n);
                        /*return vector n1 in same plane as n2,c and ca*/  
   th1=th-geo.tt[4]-geo.tt[3];                    
   planevec(nc,cca,cdn,sn,geo.dsqrt[2],th,th1,n,c);
   rotaxis(geo.de,geo.ep,-psi,cac,rom); /*rom is the rotation matrix for rotation of phi about 
       unit vector nca, vector nca must pass through origin*/
   // Note origin for psi rotation = Ca =pos[aanum][1] 
   if(ifind!=1){
     /*for cb need define a third vector to define cbeta which is out of plane*/
     for(k=0;k<=2;k++){
        cn[k]=(n[k]-pos[aanum][1][k])/geo.dsqrt[4]; /* the ca-n vector*/
     }
     crossprod(cn,cac,cc); /*cc is the component perpendicula to nca,cac*/
     for(i=0;i<=2;i++){ /* the position of cb*/
       cc[i]=cc[i]/sin(geo.tt[4]);
       cb[i]=geo.dsqrt[1]*(-0.5*cn[i]-0.5*cac[i]+sqrt(2.0/3.0)*cc[i])+pos[aanum][1][i];
     }
   /*rotate cb,n,c by psi about c-ca*/
     tranm(cb,pos[aanum][1]); /* puts Cb  at origin  = Ca*/
     multvec(rom,cb,pos[aanum][5]);
     tran(pos[aanum][5],pos[aanum][1]);
   }
   tranm(n,pos[aanum][1]);tranm(c,pos[aanum][1]);
   multvec(rom,n,pos[aanum][4]);tran(pos[aanum][4],pos[aanum][1]);
   multvec(rom,c,ct);tran(ct,pos[aanum][1]);
   /*rotate c, by phi about cn*/
   for(i=0;i<=2;i++) // find phi rotation vector  n to Ca (origin  is at n =pos[aanum][4] )
     nca[i]=(pos[aanum][1][i]-pos[aanum][4][i])/geo.dsqrt[4]; /* this is different than cn above*/
   rotaxis(geo.de,geo.ep,-phi,nca,rom);
   tranm(ct,pos[aanum][1]);
   multvec(rom,ct,pos[aanum-1][2]); 
   tran(pos[aanum-1][2],pos[aanum][1]); 
   /* find o and ca1 in same plane as c,ca2, and n*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      cca[k]=pos[aanum][4][k]-pos[aanum][1][k]; /*the ca to n vector*/
      nc[k]=pos[aanum-1][2][k]-pos[aanum][4][k]; /* the n to c vector*/
   }
   mag1=sqrt(dot(cca,cca));
   mag2=sqrt(dot(nc,nc));
   for(k=0;k<=2;k++){
      cca[k]=cca[k]/mag1;
      nc[k]=nc[k]/mag2;
   }
   cdn=dot(cca,nc);
   th=acos(cdn);/*theta= angle between 2 basis vectors n and c >0 (anticlockwise)*/
   cdn=cos(th); /* = cn dot nca = cos(theta)*/
   sn=1-cdn*cdn; /*= sqr(sin(theta)*/
   th1=th-PI+geo.tt[1];
   planevec(cca,nc,cdn,sn,geo.dsqrt[1],th,th1,pos[aanum-1][2],pos[aanum-1][1]);
   if(ifind!=1){
     th1=th+geo.tt[2]+geo.tt[1]-PI;
     planevec(cca,nc,cdn,sn,geo.dsqrt[3],th,th1,pos[aanum-1][2],pos[aanum-1][3]); 
   }     
}
#endif

void makecisrevalanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi) 
/* reverse of makealanine2: 
For cis: 2 changes geo.tt[1] = 2*PI - geo.tt[1]  and geo.tt[2] = - geo.tt[2]
This is applied to the cis pro residue itself (see makecisalanine2 for difference)    
 This is new  version - Identical results but explicitly use -psi and -phi for rotatino
   Changes:  1) Use ca to c vector  and Ca as origin for psi rotation
             2) Use n to Ca vector and Ca as origin of phi rotatin
   Use -phi and -psi because phi indicate how to form residue 2 give 1
      Here want to go from residue 2 to residue 1
   NOTE:  Any point on the axis of rotation can
     be used as origin - Thus, could use either Ca or N for
     phi rotation orign
Given  Ca1,C1 abd N2 find a.a. 0 (N1,C0,O0,Ca0, etc) if aanum=1
Does not assume ideal distances and angles in previous a.a.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2
if ifind=1, do not find o or cb
*/
{
   int i,k;
   float mag1,mag2;
   float ct[3];/*intermed co-ord*/
   float cca[3],cc[3],cn[3],nc[3],cac[3],n[3],c[3],cb[3],nca[3];
   float rom[3][3],cdn,sn,th,th1;
   float tt1,tt2;
   
//   geo.cis[aanum-1] = 1;  // Indicate is cis proline
   tt1 = 2*PI - geo.tt[1];
   tt2 = - geo.tt[2];
   /*here we set phi=psi=0 and find atoms positons*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      cca[k]=pos[aanum][1][k]-pos[aanum][2][k]; /*the c to ca vector*/
      nc[k]=pos[aanum][2][k]-pos[aanum+1][4][k]; /* the n to c vector*/
   }
   mag1=sqrt(dot(cca,cca));
   mag2=sqrt(dot(nc,nc));
   for(k=0;k<=2;k++){
      cca[k]=cca[k]/mag1;
      cac[k]= -cca[k]; /* the ca-c vector  = psi rotation vector*/
      nc[k]=nc[k]/mag2;
   }
   /*find n and c for phi=psi=0;*/
   cdn=dot(cca,nc);
   th=acos(cdn);/*theta= angle between 2 basis vectors n and c >0 (anticlockwise)*/
   sn=1-cdn*cdn; /*= sqr(sin(theta)*/
   th1=th+PI-geo.tt[4];
   planevec(nc,cca,cdn,sn,geo.dsqrt[4],th,th1,pos[aanum][1],n);
                        /*return vector n1 in same plane as n2,c and ca*/  
   th1=th-geo.tt[4]-geo.tt[3];                    
   planevec(nc,cca,cdn,sn,geo.dsqrt[2],th,th1,n,c);
   rotaxis(geo.de,geo.ep,-psi,cac,rom); /*rom is the rotation matrix for rotation of phi about 
       unit vector nca, vector nca must pass through origin*/
   // Note origin for psi rotation = Ca =pos[aanum][1] 
   if(ifind!=1){
     /*for cb need define a third vector to define cbeta which is out of plane*/
     for(k=0;k<=2;k++){
        cn[k]=(n[k]-pos[aanum][1][k])/geo.dsqrt[4]; /* the ca-n vector*/
     }
     crossprod(cn,cac,cc); /*cc is the component perpendicula to nca,cac*/
     for(i=0;i<=2;i++){ /* the position of cb*/
       cc[i]=cc[i]/sin(geo.tt[4]);
       cb[i]=geo.dsqrt[1]*(-0.5*cn[i]-0.5*cac[i]+sqrt(2.0/3.0)*cc[i])+pos[aanum][1][i];
     }
   /*rotate cb,n,c by psi about c-ca*/
     tranm(cb,pos[aanum][1]); /* puts Cb  at origin  = Ca*/
     multvec(rom,cb,pos[aanum][5]);
     tran(pos[aanum][5],pos[aanum][1]);
   }
   tranm(n,pos[aanum][1]);tranm(c,pos[aanum][1]);
   multvec(rom,n,pos[aanum][4]);tran(pos[aanum][4],pos[aanum][1]);
   multvec(rom,c,ct);tran(ct,pos[aanum][1]);
   /*rotate c, by phi about cn*/
   for(i=0;i<=2;i++) // find phi rotation vector  n to Ca (origin  is at n =pos[aanum][4] )
     nca[i]=(pos[aanum][1][i]-pos[aanum][4][i])/geo.dsqrt[4]; /* this is different than cn above*/
   rotaxis(geo.de,geo.ep,-phi,nca,rom);
   tranm(ct,pos[aanum][1]);
   multvec(rom,ct,pos[aanum-1][2]); 
   tran(pos[aanum-1][2],pos[aanum][1]); 
   /* find o and ca1 in same plane as c,ca2, and n*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      cca[k]=pos[aanum][4][k]-pos[aanum][1][k]; /*the ca to n vector*/
      nc[k]=pos[aanum-1][2][k]-pos[aanum][4][k]; /* the n to c vector*/
   }
   mag1=sqrt(dot(cca,cca));
   mag2=sqrt(dot(nc,nc));
   for(k=0;k<=2;k++){
      cca[k]=cca[k]/mag1;
      nc[k]=nc[k]/mag2;
   }
   cdn=dot(cca,nc);
   th=acos(cdn);/*theta= angle between 2 basis vectors n and c >0 (anticlockwise)*/
   cdn=cos(th); /* = cn dot nca = cos(theta)*/
   sn=1-cdn*cdn; /*= sqr(sin(theta)*/
//   th1=th-PI+geo.tt[1]; // standard trans pro
   th1 = th -PI + tt1;  // cis proline  
//printf("cis th1 = %5.3f\n",th1);
   planevec(cca,nc,cdn,sn,geo.dsqrt[1],th,th1,pos[aanum-1][2],pos[aanum-1][1]);
   if(ifind!=1){
//     th1=th+geo.tt[2]+geo.tt[1]-PI; // trans
     th1=th+tt2 + tt1 -PI;
     planevec(cca,nc,cdn,sn,geo.dsqrt[3],th,th1,pos[aanum-1][2],pos[aanum-1][3]); 
   }     
}

 
void findcb(struct geometry& geo,float n[3],float ca[3],
               float c[3],float cb[3])
/* returns the position of cb, given n,ca and c
*/
{
   int k;
   float r1,t4;
   float av[3],bv[3],crs[3],maga,magb,xx;

   r1= geo.dsqrt[1];
   t4 = geo.tt[4];
   for(k=0;k<=2;k++){
      av[k]=n[k]-ca[k];
      bv[k]=c[k]-ca[k];
   }
   maga=sqrt(dot(av,av));
   magb=sqrt(dot(bv,bv));
   for(k=0;k<=2;k++){
      av[k]=av[k]/maga;
      bv[k]=bv[k]/magb;
   }
   crossprod(av,bv,crs);
   for(k=0;k<=2;k++)
      crs[k]=crs[k]/sin(t4);
   xx=sqrt(2.0/3.0);
   for(k=0;k<=2;k++)
      cb[k]=ca[k]+r1*(-av[k]/2.0-bv[k]/2.0+xx*crs[k]);
}

  
void makealanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi) 
/* modify makealanine,
 Given C1, N2, and Ca2, find next  c2,o2,cb2,n3,ca3etc if  aanum=2
aanum is the number of last a.a, make aa=anmum+1;
NOTE: does not assume ideal length and angles
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
if ifind=1, do not find o or cb
*/
{
   int i,k;
   float uv[3],vv[3],cc[3],cb[3],can[3],cac[3],nt[3];
   float mag1,mag2,udv,sn,th1,th,c[3],n[3],rom[3][3];

//   geo.cis[aanum] = 0;  // Indicate is trans proline
   /*here we set phi=psi=0 and find positons of c and n*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      uv[k]=pos[aanum][4][k]-pos[aanum-1][2][k]; /*the C to N vector*/
      vv[k]=pos[aanum][1][k]-pos[aanum][4][k]; /* the N to Ca vector*/
   }
   mag1=sqrt(dot(uv,uv));
   mag2=sqrt(dot(vv,vv));
   for(k=0;k<=2;k++){
      uv[k]=uv[k]/mag1;
      vv[k]=vv[k]/mag2;
   }
   rotaxis(geo.de,geo.ep,phi,vv,rom); /*rot is the rotation matrix for rotation of phi about 
       unit vector nca, vector nca must pass through origin*/
   udv=dot(uv,vv); /*u dot v */
   th= -acos(udv);/*theta= angle between 2 basis vectors u and cvthis must
          correct sign so that th2=th1-th  <0 is clockwise*/
   sn=1-udv*udv; /*= sqr(sin(theta)*/
   th1=th-PI+geo.tt[4];
   planevec(uv,vv,udv,sn,geo.dsqrt[1],th,th1,pos[aanum][1],c); 
           /*returns vector c in same plane as u and v, with ca2 as origin*/
   th1=th1-PI+geo.tt[1];
   planevec(uv,vv,udv,sn,geo.dsqrt[2],th,th1,c,n); 
   if(ifind!=1){
     /*for cb need define a third vector to define cbeta which is out of plane*/
     for(i=0;i<=2;i++){
        can[i]= -vv[i]; /* the ca-n vector*/
        cac[i]=(c[i]-pos[aanum][1][i])/geo.dsqrt[1]; /* new ca-c vector*/
     }
     crossprod(can,cac,cc); /*cc is the component perpendicula to nca,cac*/
     for(i=0;i<=2;i++){ /* the geo.position of cb*/
       cc[i]=cc[i]/sin(geo.tt[4]);
       cb[i]=geo.dsqrt[1]*(-0.5*can[i]-0.5*cac[i]+sqrt(2.0/3.0)*cc[i])+pos[aanum][1][i];
     }
   /*rotate cb,c,n2 by phi about n-ca=vv with n at origin*/
     tranm(cb,pos[aanum][4]);multvec(rom,cb,pos[aanum][5]);
     tran(pos[aanum][5],pos[aanum][4]);
   } /* if ifind!=1*/
   tranm(c,pos[aanum][4]); /*translates c to c-n, puts n at origin*/
   tranm(n,pos[aanum][4]);
   multvec(rom,c,pos[aanum][2]);
   multvec(rom,n,nt);
   tran(pos[aanum][2],pos[aanum][4]);
   tran(nt,pos[aanum][4]);
   /*rotate nt by psi about cac*/
   for(i=0;i<=2;i++)
     cac[i]=(pos[aanum][2][i]-pos[aanum][1][i])/geo.dsqrt[1]; /* this is different than cac above*/
   rotaxis(geo.de,geo.ep,psi,cac,rom);
   tranm(nt,pos[aanum][1]);
   multvec(rom,nt,pos[aanum+1][4]); 
   tran(pos[aanum+1][4],pos[aanum][1]);
   /* add ca and o in same plane as ca-c and c-n*/
   for(k=0;k<=2;k++){
     uv[k]=(pos[aanum][2][k]-pos[aanum][1][k])/geo.dsqrt[1];
     vv[k]=(pos[aanum+1][4][k]-pos[aanum][2][k])/geo.dsqrt[2];
   }
   th= -PI+geo.tt[1];
   th1=th+PI-geo.tt[3];
   udv=cos(th); /* = cn dot nca = cos(theta)*/
   sn=1-udv*udv; /*= sqr(sin(theta)*/
   planevec(uv,vv,udv,sn,geo.dsqrt[4],th,th1,pos[aanum+1][4],pos[aanum+1][1]); 
   if(ifind!=1){
      th1=PI-geo.tt[2];   
      planevec(uv,vv,udv,sn,geo.dsqrt[3],th,th1,pos[aanum][2],pos[aanum][3]);
   }    
}

  
void makecisalanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi) 
/*
NOTE: This must be used for residue preceding Cis proline because
   the last Ca is the only thing that is changed and this is the
   the Ca in the next residue.
Used for Cis proline  
  ONLY change geo.tt[3] (angle between C - N and new N-Ca) to
    2PI - geo.tt[3]
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
if ifind=1, do not find o or cb
*/
{
   int i,k;
   float uv[3],vv[3],cc[3],cb[3],can[3],cac[3],nt[3];
   float mag1,mag2,udv,sn,th1,th,c[3],n[3],rom[3][3];

   geo.cis[aanum] = 1;  // Indicate is cis proline
   /*here we set phi=psi=0 and find positons of c and n*/
   for(k=0;k<=2;k++) { /*these are the known initial (unit) vectors*/
      uv[k]=pos[aanum][4][k]-pos[aanum-1][2][k]; /*the C to N vector*/
      vv[k]=pos[aanum][1][k]-pos[aanum][4][k]; /* the N to Ca vector*/
   }
   mag1=sqrt(dot(uv,uv));
   mag2=sqrt(dot(vv,vv));
   for(k=0;k<=2;k++){
      uv[k]=uv[k]/mag1;
      vv[k]=vv[k]/mag2;
   }
   rotaxis(geo.de,geo.ep,phi,vv,rom); /*rot is the rotation matrix for rotation of phi about 
       unit vector nca, vector nca must pass through origin*/
   udv=dot(uv,vv); /*u dot v */
   th= -acos(udv);/*theta= angle between 2 basis vectors u and cvthis must
          correct sign so that th2=th1-th  <0 is clockwise*/
   sn=1-udv*udv; /*= sqr(sin(theta)*/
   th1=th-PI+geo.tt[4];
   planevec(uv,vv,udv,sn,geo.dsqrt[1],th,th1,pos[aanum][1],c); 
           /*returns vector c in same plane as u and v, with ca2 as origin*/
   th1=th1-PI+geo.tt[1];
   planevec(uv,vv,udv,sn,geo.dsqrt[2],th,th1,c,n); 
   if(ifind!=1){
     /*for cb need define a third vector to define cbeta which is out of plane*/
     for(i=0;i<=2;i++){
        can[i]= -vv[i]; /* the ca-n vector*/
        cac[i]=(c[i]-pos[aanum][1][i])/geo.dsqrt[1]; /* new ca-c vector*/
     }
     crossprod(can,cac,cc); /*cc is the component perpendicula to nca,cac*/
     for(i=0;i<=2;i++){ /* the geo.position of cb*/
       cc[i]=cc[i]/sin(geo.tt[4]);
       cb[i]=geo.dsqrt[1]*(-0.5*can[i]-0.5*cac[i]+sqrt(2.0/3.0)*cc[i])+pos[aanum][1][i];
     }
   /*rotate cb,c,n2 by phi about n-ca=vv with n at origin*/
     tranm(cb,pos[aanum][4]);multvec(rom,cb,pos[aanum][5]);
     tran(pos[aanum][5],pos[aanum][4]);
   } /* if ifind!=1*/
   tranm(c,pos[aanum][4]); /*translates c to c-n, puts n at origin*/
   tranm(n,pos[aanum][4]);
   multvec(rom,c,pos[aanum][2]);
   multvec(rom,n,nt);
   tran(pos[aanum][2],pos[aanum][4]);
   tran(nt,pos[aanum][4]);
   /*rotate nt by psi about cac*/
   for(i=0;i<=2;i++)
     cac[i]=(pos[aanum][2][i]-pos[aanum][1][i])/geo.dsqrt[1]; /* this is different than cac above*/
   rotaxis(geo.de,geo.ep,psi,cac,rom);
   tranm(nt,pos[aanum][1]);
   multvec(rom,nt,pos[aanum+1][4]); 
   tran(pos[aanum+1][4],pos[aanum][1]);
   /* add ca and o in same plane as ca-c and c-n*/
   for(k=0;k<=2;k++){
     uv[k]=(pos[aanum][2][k]-pos[aanum][1][k])/geo.dsqrt[1];
     vv[k]=(pos[aanum+1][4][k]-pos[aanum][2][k])/geo.dsqrt[2];
   }
   th= -PI+geo.tt[1];  // standard
//   th1=th+PI-geo.tt[3]; // standard
   th1 = th -PI+ geo.tt[3];  // cis
   udv=cos(th); /* = cn dot nca = cos(theta)*/
   sn=1-udv*udv; /*= sqr(sin(theta)*/
   planevec(uv,vv,udv,sn,geo.dsqrt[4],th,th1,pos[aanum+1][4],pos[aanum+1][1]); 
   if(ifind!=1){
      th1=PI-geo.tt[2];   
      planevec(uv,vv,udv,sn,geo.dsqrt[3],th,th1,pos[aanum][2],pos[aanum][3]);
   }    
}


void addcbn(struct geometry& geo,float pos[NAFIT][NAA][3],int ifit,int iatom)
/*
 Add the i=iatom (>=3) side chain atom to the ifit residue
 Procedure:
   Similar to addcb1:
    Add cb[i] in the cb[i-3] -cb[i-2]  - cv[i-1] plane
  If iatom = 3 then use Ca instead of cb[i-3]
*/
{
   int k,i1,i2,i3;
   float nv[3],cv[3],rdn,ndc,det,rdc,be,al,rv[3],mag1,mag2; 

   if(iatom ==3)
      i1 = 1; // Ca
   else
      i1 = iatom -3 + 4;
   i2 = iatom-2+4;
   i3 = iatom-1+4;  
   for(k=0;k<=2;k++){
      nv[k]=(pos[ifit][i2][k]-pos[ifit][i1][k]); /*cb[i-3] to cb[i-2] vector*/
      cv[k]=(pos[ifit][i3][k]-pos[ifit][i2][k]); /*cb[i-2] to cb[i-1] vector*/
   }
   mag1 = sqrt(dot(nv,nv));// mag1 and mag2 should = dsqrt[4]
   mag2 = sqrt(dot(cv,cv));
   for(k=0;k<=2;k++){
      nv[k] = nv[k]/mag1;
      cv[k] = cv[k]/mag2;
   }
   rdn=1.0; /* r dot n*/
   rdc=ndc=cos(PI-geo.tt[4]);/* r dot c and n dot c */
   det=1-ndc*ndc;
   al=(rdn-rdc*ndc)/det;
   be=(rdc-rdn*ndc)/det;
   for(k=0;k<=2;k++){
     rv[k]=al*nv[k]+be*cv[k];
     pos[ifit][4+iatom][k]=pos[ifit][4+iatom-1][k]+geo.dsqrt[4]*rv[k]; /*the second sidechain atom*/
   } 
}
   

void addcb2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifit)
 /* Add second sidechain atom in N-Ca-Cb plane: geo.pos[6]=cb+r4*rv[]*/
{
   int k;
   float nv[3],cv[3],rdn,ndc,det,rdc,be,al,rv[3],mag1,mag2;
   for(k=0;k<=2;k++){
      nv[k]=(pos[ifit][1][k]-pos[ifit][4][k]); /*n to Ca vector*/
      cv[k]=(pos[ifit][5][k]-pos[ifit][1][k]); /*Ca to Cb vector*/
   }
   mag1 = sqrt(dot(nv,nv));// mag1 and mag2 should = dsqrt[4]
   mag2 = sqrt(dot(cv,cv));
   for(k=0;k<=2;k++){
      nv[k] = nv[k]/mag1;
      cv[k] = cv[k]/mag2;
   }
   rdn=1.0; /* r dot n*/
   rdc=ndc=cos(PI-geo.tt[4]);/* r dot c and n dot c */
   det=1-ndc*ndc;
   al=(rdn-rdc*ndc)/det;
   be=(rdc-rdn*ndc)/det;
   for(k=0;k<=2;k++){
     rv[k]=al*nv[k]+be*cv[k];
     pos[ifit][6][k]=pos[ifit][5][k]+geo.dsqrt[4]*rv[k]; /*the second sidechain atom*/
   } 
}

void maketestpep(struct geometry& geo,float pos[NAFIT][NAA][3],int numaa,
        float phi,float psi)
/* returns test peptide of numaa amino acids,  beginning at geo.stpos
    all  phi psi (in radians)  are the same 
phi/psi[0,1,2,  ] are the array of angles
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
For helix:
  phi= -60.0*DEGTORAD;  value used in for helix 
   psi= -45.0*DEGTORAD;
 and rev angles: phirev = 135, psirev = -120
For sheet:

*/
{
   int i,k;
   firstaa2(geo); /*returns a.a (caw, etc) with ca1 at origin*/
   for(k=0;k<=2;k++){ 
      pos[1][1][k]=geo.caw[1][k]+geo.stpos[k];
      pos[1][2][k]=geo.cw[1][k]+geo.stpos[k];
      pos[1][3][k]=geo.ow[1][k]+geo.stpos[k];
      pos[2][4][k]=geo.nw[2][k]+geo.stpos[k];
      pos[2][1][k]=geo.caw[2][k]+geo.stpos[k];
   }
   /*at this point should have one a.a (ca1 to ca2) using geo.pos variables*/
   for(i=2;i<=numaa;i++) 
      makealanine2(geo,pos,0,i,phi,psi);
//   geo.nend[0]=1; // not used now
//   geo.nend[1]=numaa+1;
   geo.nafit=numaa+1;
   if(geo.nside==2){ // add cb2
     for(i=2;i<geo.nafit;i++){ /* sidechain not defined at two ends*/
        addcb2(geo,pos,i); /*add second carbon on sidechain*/
        geo.iside[i]=2;
     }
   }
   else{ // geo.nside =1
     for(i=2;i<geo.nafit;i++) /* sidechain not defined at two ends*/
        geo.iside[i]=1;
   }     
   /*at this point endnum[0] and endnum[1] are are close to correct values
     these are checked after every time step*/
} 

void findnextpos(struct geometry& geo,int resi,int ai2,
    int *resf,int *af2,int st[3],float len,float posi[3],float posf[3])
/*  Finds the next position at dist len from intial posi
     posi and posf are initial and final pos
     initial pos is less then atom ai2 on resi 
     final pos in less then atom af2 on resf)
Only look forward distance of two atoms - if greater than that - get error and exit routine
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i][res][x,y,z]...
  define resi = 0 (Ca), 1 (C) or 2 (N next residue)
  and convert as follows:
    st[0]=0;  Ca positon  Forward direction : 0 to 1 to 2  (Ca to C to N of next residue)
    st[1]=1;  C 
    st[2] = 4;  N
If need to go to next atom  - use quadrtic equation:
  vec2 = vector from atom ai2 to next atom
  vec = vector from posi to ai2
  xter = fractional distance along vec2 that give posf
  pos = vec + x *vec2
  quadratic equaton:  x*x (dot(vec2,vec2)  + x (2*dot(vec,vec2) + dot(vec,vec) = len*len  
*/
{
    int i,nextres,nextat;
    float vec[3],vec2[3],vec3[3],dist2,dist1,xter;
    float aa,bb,cc;

    for(i= 0;i<=2;i++)
       vec[i] = geo.pos[resi][ st[ai2] ][i] - posi[i];
    dist1 = sqrt( dot(vec,vec) ); // dist from posi to next atom
    if(len < dist1){
       for(i=0;i<=2;i++)
         posf[i] = posi[i]+vec[i]*len/dist1;
       *af2 = ai2;
       *resf= resi;
    }
    else{ // len>=dist1, find vec to next atom
      if( ai2 == 1)
        nextres = resi+1;
      else
        nextres = resi;
      nextat = (ai2+1) % 3;  // Next atom modulo e
      for(i= 0;i<=2;i++) // vector from ai2 to nextat
        vec2[i] = geo.pos[nextres][ st[nextat] ][i] - geo.pos[resi][ st[ai2] ][i];
      for(i= 0;i<=2;i++) 
        vec3[i] = geo.pos[nextres][ st[nextat] ][i] - posi[i];
      dist2 = sqrt( dot(vec3,vec3) ); // dist from next atom to initial pos
      if(len < dist2){ // Then solve quadratic to find distance xter along vec that equal len
        cc= sqr(dist1)-sqr(len); // should be less than 0
        bb = 2*dot(vec,vec2);
        aa = dot(vec2,vec2);
        xter = (-bb+sqrt(bb*bb -4.0*aa*cc) )/(2.0*aa); // = fractional distance along bond towrard nextatom
        for(i=0;i<=2;i++)
          posf[i] =posi[i]+ vec[i]+xter*vec2[i];
        *af2 = nextat;
        *resf = nextres;
      }
      else{ //len>= dist1  - must goto next atom (total of 2 atoms)
        ai2= nextat;
        dist1=dist2;
        resi = nextres;
        for(i=0;i<=2;i++)
            vec[i] = vec2[i];
        if( ai2 == 1)
          nextres = resi+1;
         else
           nextres = resi;
         nextat = (nextat+1) % 3;  // Next atom modulo e
         for(i= 0;i<=2;i++) // vector from ai2 to nextat
           vec2[i] = geo.pos[nextres][ st[nextat] ][i] - geo.pos[resi][ st[ai2] ][i];
         for(i= 0;i<=2;i++) 
           vec3[i] = geo.pos[nextres][ st[nextat] ][i] - posi[i];
         dist2 = sqrt( dot(vec3,vec3) ); // dist from next atom to initial pos
         if(len < dist2){ // Then solve quadratic to find distance xter along vec that equal len
           aa= sqr(dist1)-sqr(len); // should be less than 0
           bb = 2*dot(vec,vec2);
           cc = dot(vec2,vec2);
           xter = (-bb+sqrt(bb*bb -4.0*aa*cc) )/(2.0*aa); // = fractional distance along bond towrard nextatom
           for(i=0;i<=2;i++)
             posf[i] =posi[i]+ vec[i]+xter*vec2[i];
           *af2 = nextat;
           *resf = nextres;
         }
         else{
           fprintf(fpout,"ERROR: Too big a step in findnextpos - indicates grid size is very large\n");
           exit(0);
         }
      } // second else
   } // first else
}
        
void findprevpos(struct geometry& geo,int resi,int ai2,
    int *resf,int *af2,int st[3],float len,float posi[3],float posf[3])
/*  Finds the previous position at dist len from intial posi
     posi and posf are initial and final pos
     initial pos is greater than atom ai2 on resi 
     final pos in greater then atom af2 on resf)
Only look backward distance of two atoms - if greater than that - get error and exit routine
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i][res][x,y,z]...
  define resi = 0 (Ca), 1 (C) or 2 (N next residue)
  and convert as follows:
    st[0]=0;  Ca positon  Forward direction : 0 to 1 to 2  (Ca to C to N of next residue)
    st[1]=1;  C 
    st[2] = 4;  N
If need to go to next atom  - use quadrtic equation:
  vec2 = vector from atom ai2 to previous atom
  vec = vector from posi to ai2
  xter = fractional distance along vec2 that give posf
  pos = vec + x *vec2
  quadratic equaton:  x*x (dot(vec2,vec2)  + x (2*dot(vec,vec2) + dot(vec,vec) = len*len  
*/
{
    int i,prevres,prevat;
    float vec[3],vec2[3],vec3[3],dist2,dist1,xter;
    float aa,bb,cc;

    for(i= 0;i<=2;i++)
       vec[i] = geo.pos[resi][ st[ai2] ][i] - posi[i];
    dist1 = sqrt( dot(vec,vec) ); // dist from posi to next atom
    if(len < dist1){
       for(i=0;i<=2;i++)
         posf[i] = posi[i]+vec[i]*len/dist1;
       *af2 = ai2;
       *resf= resi;
    }
    else{ // len>=dist1, find vec to next atom
      if( ai2 == 2)  // N
        prevres = resi-1;
      else
        prevres = resi;
      prevat = ai2-1;
      if(prevat < 0)
        prevat = prevat+3;  // This is equivalent to modulo procedure
      for(i= 0;i<=2;i++) // vector from ai2 to prevat
        vec2[i] = geo.pos[prevres][ st[prevat] ][i] - geo.pos[resi][ st[ai2] ][i];
      for(i= 0;i<=2;i++) 
        vec3[i] = geo.pos[prevres][ st[prevat] ][i] - posi[i];
      dist2 = sqrt( dot(vec3,vec3) ); // dist from next atom to initial pos
      if(len < dist2){ // Then solve quadratic to find distance xter along vec that equal len
        cc= sqr(dist1)-sqr(len); // should be less than 0
        bb = 2*dot(vec,vec2);
        aa = dot(vec2,vec2);
        xter = (-bb+sqrt(bb*bb -4.0*aa*cc) )/(2.0*aa); // = fractional distance along bond towrard nextatom
        for(i=0;i<=2;i++)
          posf[i] =posi[i]+ vec[i]+xter*vec2[i];
        *af2 = prevat;
        *resf = prevres;
      }
      else{ //len>= dist1  - must goto next atom (total of 2 atoms)
        ai2= prevat;
        dist1=dist2;
        resi = prevres;
        for(i=0;i<=2;i++)
            vec[i] = vec2[i];
        if( ai2 == 2)
          prevres = resi-1;
         else
           prevres = resi;
         prevat = prevat-1;
         if(prevat < 0)
           prevat = prevat+3;  // This is equivalent to modulo procedure
         for(i= 0;i<=2;i++) // vector from ai2 to nextat
           vec2[i] = geo.pos[prevres][ st[prevat] ][i] - geo.pos[resi][ st[ai2] ][i];
         for(i= 0;i<=2;i++) 
           vec3[i] = geo.pos[prevres][ st[prevat] ][i] - posi[i];
         dist2 = sqrt( dot(vec3,vec3) ); // dist from next atom to initial pos
         if(len < dist2){ // Then solve quadratic to find distance xter along vec that equal len
           aa= sqr(dist1)-sqr(len); // should be less than 0
           bb = 2*dot(vec,vec2);
           cc = dot(vec2,vec2);
           xter = (-bb+sqrt(bb*bb -4.0*aa*cc) )/(2.0*aa); // = fractional distance along bond towrard nextatom
           for(i=0;i<=2;i++)
             posf[i] =posi[i]+ vec[i]+xter*vec2[i];
           *af2 = prevat;
           *resf = prevres;
         }
         else{
           fprintf(fpout,"ERROR: Too big a step in findprevpos - indicates grid size is very large\n");
           exit(0);
         }
      } // second else
   } // first else
}

#if 0
void findallowedang(struct griddata& den1,struct geometry& geo)
/*
THIS IS NOT USED NOW - INSTEAD USE ROUTIEN TO ESTIMATE PHI PSI 'estphipsi'
Finds allowed range of theta and phi for fourt mainch point, given 3 pts
   Define av = unit vector from first to second pt.
          bv    ..              second to third
          cv    ....            thrid to fourt
Old   Set up coordinate system with z = bv; x = z cross av (make unit lenght) and y = z cross ax
Since Ca - C - N -Ca all lie in same plane -
   I start at Ca or C, all 3 vectors will be in same plane, thus seems
      natural to define z as perpendicular to this plane, so that
      theta = 0;  Thus:
New   Set up coordinate system with x = av; z = av cross bv (make unit lenght)
         and y = zv cross av

   Then write cv in theta, phi coord:  cos(theta) = cv dot z, and
      tan(phi) = (cv dot y)/(cv dot x);
  Because of planar nature of bond, most theta angle should = PI
    Thus consider two cases:
    I) if theta = 0,  phimin,phimax   forhelang[p1][p2][p3][0,1]
    II) if theta not PI, then 
          phimin,phimax   given by forhelang[p1][p2][p3][2,3]
          thetamin,thetamax given by forhelang[p1][p2][p3][4,5]
Max and min are found by starting with initial point at Ca, C and N position
  and intermediate postiitions respectivel
The possible value of len1 adn len2 are determined by
   grid dimension and type of mainch step (i.e.   i,j,k  to i+1,j,k   or i,j,k to i+1,j+1,k+1  etc
There are a total of  7  lengths:
  characterized by paramter  p = delta i  + 2* delta j + 4* delta k -1   (where delta is absolute value)
                           p value
   i,j,k   to i+1,j,k         0    (Note, i,j,k to i-1,j,k is identical to this)
   i,j,k      i,j+1,k         1
   i,j,k      i,j,k+1         3
   i,j,k      i+1,j+1,k       2
   i,j,k      i+1,j, k+1      4
   i,j,k      i,j+1,k+1       5
   i,j,k      i+1,j+1,k+1     6
Procedure:  Starting at Ca, C or N and half way posstions
    run through all lengts and calculate theta and phi - find range
  Do this for both forward and backward direction
dsqrt[1]= Ca - C distance; dsqrt[2] = C - N; dsqrt[4] = Ca - N distance
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
*/
{
   int i,ii,j,k,p,ivec[3],st[3],p1,p2,p3;
   int resi1,resi2,resf,ai21,ai22,af2;
   int sres[10],satom[10];
   float spos[10][3],ca2[3];
   float len[8],posf[3],posi1[3],posi2[3],phi,theta;
   float dist,av[3],bv[3],cv[3],xx[3],yy[3],zz[3],fvec[3],fvec0[3];
   float del = PI/2000.0;

   fprintf(fpout," Finding allowed angle range for all posible steps len1,len2, and len3\n");
   ivec[0]=ivec[1]=ivec[2]=0;
   imcv(den1,ivec,fvec0);// convert int grid ivec to float world co-ord fvec
   for(i=0;i<=1;i++) // find lengths as function of p
     for(j=0;j<=1;j++)
       for(k=0;k<=1;k++){
         p=i +2*j +4*k -1;
         if(p>=0){
            ivec[0]=i;ivec[1]=j;ivec[2]=k;
            imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
            len[p] = distv(fvec0,fvec);
         }
       }
    st[0]=1; // Ca positon  Forward direction : 0 to 1 to 2  (Ca to C to N of next residue)
    st[1]=2; // C 
    st[2] = 4; // N
    for(p1=0;p1<=6;p1++) // set initial phi and theta range
      for(p2=0;p2<=6;p2++)
        for(p3=0;p3<=6;p3++){
           geo.forhelang[p1][p2][p3][0] = 1000.0; //set minimum to big value
           geo.forhelang[p1][p2][p3][1]= -1000.0; //set maximum to small value
           geo.forhelang[p1][p2][p3][2]= 1000.0; //set minimum to big value
           geo.forhelang[p1][p2][p3][3] = -1000.0; //set maximum to small value
           geo.forhelang[p1][p2][p3][4] = 1000.0; //set minimum to big value
           geo.forhelang[p1][p2][p3][5] = -1000.0; //set maximum to small value
        }
// Set up test start positons at Ca, C and N and intermediate positons
    for(i=0;i<=2;i++){
      spos[0][i] = geo.pos[1 ][1][i]; // Ca
      spos[2][i] = geo.pos[1][2 ][i]; // C
      spos[4][i] = geo.pos[2][4 ][i]; // N
      spos[1][i] = (spos[0][i]+spos[2][i])/2.0;
      spos[3][i] = (spos[2][i]+spos[4][i])/2.0;
      ca2[i] = geo.pos[2][1 ][i]; // Next  Ca
      spos[5][i] = (spos[4][i]+ca2[i])/2.0;
    }
    sres[0]=sres[1]=1; // The residue number of the next atom position
    sres[2]=sres[3]=sres[4]=sres[5]=2;
    satom[0]=satom[1] = 1; // The next atom number  (0,1,or 2)
    satom[2]=satom[3] = 2;
    satom[4] = satom[5]= 0;
    for( ii = 0;ii<=5;ii++){ // run though different starting atom positions
       for(p1=0;p1<=6;p1++){
         findnextpos(geo,sres[ii],satom[ii],&resf,&af2,st,len[p1],spos[ii],posf);
         for(i=0;i<=2;i++)
            av[i]= posf[i]-spos[ii][i];
         for(i=0;i<=2;i++) // make unit av vector
            av[i] = av[i]/len[p1]; // len[p1] = length of av

         for(i=0;i<=2;i++) // make unit av vector
            posi1[i] = posf[i];
         resi1 = resf;
         ai21 = af2;
         for(p2=0;p2<=6;p2++){
           findnextpos(geo,resi1,ai21,&resf,&af2,st,len[p2],posi1,posf);
           for(i=0;i<=2;i++)
              bv[i]= posf[i]-posi1[i];
//New   Set up coordinate system with x = av; z = av cross bv (make unit lenght)
//         and y = zv cross av
           for(i=0;i<=2;i++) 
              xx[i] = av[i]; // already a unit vector
           crossprod(xx,bv,zz);  // make  z vector 
           crossprod(zz,xx,yy);  // make  y vector 
           dist = sqrt( dot(zz,zz) );
           for(i=0;i<=2;i++) 
              zz[i] = zz[i]/dist;
           dist = sqrt( dot(yy,yy) );
           for(i=0;i<=2;i++) 
              yy[i] = yy[i]/dist;
           for(i=0;i<=2;i++) // reset posi
              posi2[i] = posf[i];
           resi2 = resf;
           ai22 = af2;
           for(p3=0;p3<=6;p3++){
              findnextpos(geo,resi2,ai22,&resf,&af2,st,len[p3],posi2,posf);
              for(i=0;i<=2;i++)
                cv[i]= posf[i]-posi2[i];
              for(i=0;i<=2;i++) // make unit cv vector
                cv[i] = cv[i]/len[p3];

              theta = acos( dot(cv,zz) ); // range 0 to PI
              phi = atan2( dot(cv,yy),dot(cv,xx) );//atan(y,x) return arctan(y/x)  range -PI to PI
if( (p1==2)&&(p2==0)&&(p3==4) )
  fprintf(fpout,"");
              if( (theta > PI/2.0 -del)&&(theta < PI/2.0 +del) ){ // theta = PI/2.0;
                  if(phi < geo.forhelang[p1][p2][p3][0])
                      geo.forhelang[p1][p2][p3][0] = phi;
                  if(phi> geo.forhelang[p1][p2][p3][1])
                      geo.forhelang[p1][p2][p3][1] = phi;
              }
              else{ // theta not equal PI/2
                  if(phi < geo.forhelang[p1][p2][p3][2])
                      geo.forhelang[p1][p2][p3][2] = phi;
                  if(phi> geo.forhelang[p1][p2][p3][3])
                      geo.forhelang[p1][p2][p3][3] = phi;
                  if(theta < geo.forhelang[p1][p2][p3][4])
                     geo.forhelang[p1][p2][p3][4] = theta;
                  if(theta> geo.forhelang[p1][p2][p3][5])
                     geo.forhelang[p1][p2][p3][5] = theta;
             }
           }//p3
         } // p2
      } //p1
   }// ii  (starting atom)
   fprintf(fpout," Done\n");
#if 1
    for(p1=0;p1<=6;p1=p1+2) // set initial phi and theta range
      for(p2=0;p2<=6;p2=p2+2)
        for(p3=0;p3<=6;p3=p3+2){
          fprintf(fpout,"forhelang[%d][%d][%d][0 to 5] = %6.3f %6.3f  %6.3f %6.3f  %6.3f  %6.3f\n",
                 p1,p2,p3,geo.forhelang[p1][p2][p3][0],geo.forhelang[p1][p2][p3][1],
                 geo.forhelang[p1][p2][p3][2],geo.forhelang[p1][p2][p3][3],
                 geo.forhelang[p1][p2][p3][4],geo.forhelang[p1][p2][p3][5]);
        }
#endif
}
#endif


float b2funcn(float x,head2n)
/*
NEW VERSION:  arbitrary independent variable = ix, x = cc[ix]
  this is the 4th Non-linear equation
  find root by zbrent, 
  function of form  b2func(x = b2, head) = 0.0

*/
{
    int i,j;   
    float d[3],lenw1,lenw2;
    float value;  // want make value =0

    cc[ix] = x;
    for(i=0;i<=2;i++)
       d[i] = -x*r[ix][i]+dp[i];
   for(j=0;j<=2;j++){ // define each column  m[i][j]
     for(i= 0;i<=2;i++)
       m[j][i] = r[ t[i] ][j];
   }
    for(j=0;j<=2;j++) // determinant for cc[ t[0] ]
       m[j][0] = d[j];  // change first column
    cc[ t[0] ] = determ(m)/det;
    for(j=0;j<=2;j++){ // determinant for cc[ t[1] ]
       m[j][0] = r[ t[0] ][j]; // set first column back to original values
       m[j][1] = d[j];  // change second column
    }
    cc[ t[1] ] = determ(m)/det;
    for(j=0;j<=2;j++){ // determinant for cc[ t[2] ]
       m[j][1] = r[ t[1] ][j]; // set second column back to original values
       m[j][2] = d[j];  // change second column
    }
    cc[ t[2] ] = determ(m)/det;
    // Now use these values in non-linear equation with x = bb[1] = b2
// NOTE: because r[2] = -u2 and r[3] = -v2  must change sign of  some dot product
    lenw1 = sqrt( cc[0]*cc[0]*u1du1 + 2.0*cc[0]*cc[1]*u1dv1+cc[1]*cc[1]*v1dv1);
    lenw2 = sqrt( cc[2]*cc[2]*u2du2 + 2.0*cc[2]*cc[3]*u2dv2+cc[3]*cc[3]*v2dv2);
//fprintf(fpout," lenw1 = %f  lenw2 = %f\n",lenw1,lenw2);
    value = cc[0]*cc[2]*u1du2+cc[3]*cc[0]*u1dv2+cc[1]*cc[2]*v1du2+ cc[3]*cc[1]*v1dv2 -
          lenw1*lenw2*cosang;
    return value;
}


float b2func(float x,head2)
/*  this is the 4th Non-linear equation
  find root by zbrent, 
  function of form  b2func(x = b2, head) = 0.0
*/
{
    int i,j;   
    float d[3],lenw1,lenw2;
    float value;  // want make value =0

    for(i=0;i<=2;i++)
       d[i] = x*v2[i]+dp[i];
#if 1  // must redefine m each time
    for(j=0;j<=2;j++){ // define each column  m[i][j]
       m[j][0] = u1[j];
       m[j][1] = v1[j];
       m[j][2] = -u2[j];
    }     
//    det = determ(m);
//fprintf(fpout," in b2fjnc det = %f\n",det);
#endif
//fprintf(fpout,"x= %f det = %f\n",x,det);
    for(j=0;j<=2;j++) // determinant for a1
       m[j][0] = d[j];  // change first column
    aa[0] = determ(m)/det;
    for(j=0;j<=2;j++){ // determinant for b1
       m[j][0] = u1[j]; // set first column back to original values
       m[j][1] = d[j];  // change second column
    }
    bb[0] = determ(m)/det;
    for(j=0;j<=2;j++){ // determinant for a2
       m[j][1] = v1[j]; // set second column back to original values
       m[j][2] = d[j];  // change second column
    }
    aa[1] = determ(m)/det;
    // Now use these values in non-linear equation with x = bb[1] = b2
    lenw1 = sqrt( aa[0]*aa[0]*u1du1 + 2.0*aa[0]*bb[0]*u1dv1+bb[0]*bb[0]*v1dv1);
    lenw2 = sqrt( aa[1]*aa[1]*u2du2 + 2.0*aa[1]*x*u2dv2+x*x*v2dv2);
//fprintf(fpout," lenw1 = %f  lenw2 = %f\n",lenw1,lenw2);
    value = aa[0]*aa[1]*u1du2+x*aa[0]*u1dv2+bb[0]*aa[1]*v1du2+ x*bb[0]*v1dv2 -
          lenw1*lenw2*cosang;
    bb[1] = x;
    return value;
}
    

#define ITMAX 100
#define EPS 3.0e-8

float zbrent(float (*func)(float x,head2),float x1,
        float x2,head2,float tol)
{
	int iter;
	float a=x1,b=x2,c,d,e,min1,min2;
	float fa=(*func)(a,func2),fb=(*func)(b,func2),fc,p,q,r,s,tol1,xm;
//        void nrerror(char error_text[]);

	if (fb*fa > 0.0) nrerror("Root must be bracketed in ZBRENT");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (fabs(fc) < fabs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*fabs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (fabs(xm) <= tol1 || fb == 0.0) return b;
		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=fabs(p);
			min1=3.0*xm*q-fabs(tol1*q);
			min2=fabs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (fabs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
		fb=(*func)(b,func2);
	}
	nrerror("Maximum number of iterations exceeded in ZBRENT");
        return -1000.0;  // indicates error
}
  
#if 1 
void cajunct(float m[3][3],float det,float ang,float u1[3],float v1[3],float p1[3],
     float u2[3],float v2[3],float p2[3],float w1[3],float w2[3])
/* routine to find the Ca position connnects two adjacent
  peptide planes.
Given vectors  u1 and v1 in plane 1 and u2 , v2 in plane 2
  v1 and v2 end at point p1,p2.
  find  w1 and w2  such that p1+w1 = Ca = p2 + w2
Using fact that w1 must be in u1,v1 plane: w1 = a u1 + b v1
   and w2 in u2,v2 plane and angle (w1 dot w2) must = ang;
 Can uniquely find Ca.
Procedure:  4 varialbles,  a1,b1,a2,b2,
   Have 3 linear equations and 1 non-linear
  Solve the 3 linear equaatons for a1,b1,a2 as function of b2
  then Insert these values into 4th eq. and find root (x = b2)
    using zbrent. (should be bounded by 0< b2 < 2)
*/
{
   int i;
   float aa[2],bb[2],cosang,b2;
   float  u1du1,u1dv1,v1dv1;  // u1 dot u1 etc.
   float  u2du2,u2dv2,v2dv2;  // u1 dot u1 etc.
   float u1du2,u1dv2,v1du2,v1dv2,dp[3],tol;

   tol=1.0e-4;
   u2du2 = dot(u2,u2);
   u2dv2 = dot(u2,v2);
   v2dv2 = dot(v2,v2);
   u1du1 = dot(u1,u1);
   u1dv1 = dot(u1,v1);
   v1dv1 = dot(v1,v1);
   u1du2 = dot(u1,u2);
   u1dv2 = dot(u1,v2);
   v1du2 = dot(v1,u2);
   v1dv2 = dot(v1,v2);
   for(i = 0;i<=2;i++)
     dp[i] = p2[i] - p1[i];
   cosang = cos(ang);
#if 1
   for(i= 0;i<=100;i++){
      float xi = (i-50.0)/25.0;
      fprintf(fpout,"b2func( %f) = %f\n",xi,b2func(xi,func2) );
   }
#endif
   // Find bounding points - Know answer must be  >0 and < 2, start at 0, step 0.1

   float x0 = 0.0;
   float xstep = 0.05;
   float y0 = b2func(x0,func2);
   float y1 = b2func(x0+xstep,func2);
fprintf(fpout,"y0 = %f  y1 = %f\n",y0,y1);
   int isign,iplus;
   float xlow,xhigh;
   if(y0<y1){
       if(y0<0){ // increase  x1 in steps till sign changes
         isign =1;
         iplus= -1;
       }
       else{ // decrease  x1 in steps till sign changes
         isign = -1;
         iplus=1;
       }
   }
   else{ // y0>y1
       if(y0<0){ // decrease  x1 in steps till sign changes
         isign = -1;
         iplus= -1;
       }
       else{ // increase  x1 in steps till sign changes
         isign = +1;
         iplus=1;
       }
   }
   float x1 = x0;
fprintf(fpout,"isign = %d  iplus = %d\n",isign,iplus);
   do{// when exit this loop, have changed sign
     x1 = x1+isign*xstep;
     y1 = b2func(x1,func2);
     if(x1>2.0){
       fprintf(fpout,"ERROR:  CANNOT FIND 0 FOR zbrent in cajunct\n");
       exit(0);
     }
   }while(iplus*y1>0); 
   if(isign>0){
       xlow = x0;
       xhigh = x1;
   }
   else{
       xlow = x1; 
       xhigh = x0;
   }
   b2=zbrent(b2func,xlow,xhigh,func2,tol);
fprintf(fpout," At end of zbrent: a[0]= %f a[1] = %f b[0] = %f b[1] = %f  value = %f\n",
       aa[0],aa[1],bb[0],bb[1],b2func(b2,func2) );
   for(i=0;i<=2;i++){
      w1[i] = aa[0]*u1[i]+bb[0]*v1[i];
      w2[i] =  aa[1]*u2[i]+bb[1]*v2[i];
   }
#if 1 // check answer
float chpos[3];
   for(i=0;i<=2;i++)
      chpos[i] = p1[i]+w1[i];
   fprintf(fpout,"Ca theory (w1) pos =  %f %f %f  \n",chpos[0],chpos[1],chpos[2]);
   for(i=0;i<=2;i++)
      chpos[i] = p2[i]+w2[i];
   fprintf(fpout,"Ca theory (w2) pos =  %f %f %f  \n",chpos[0],chpos[1],chpos[2]);
   fprintf(fpout,"   aa[0] = %f aa[1]= %f  bb[0] = %f bb[1] = %f\n",aa[0],aa[1],bb[0],bb[1]);
#endif
}
#endif   


float zbrentn(float (*func)(float x,head2n),float x1,
        float x2,head2n,float tol)
// New version - used with cajunctn
{
	int iter;
	float a=x1,b=x2,c,d,e,min1,min2;
	float fa=(*func)(a,func2n),fb=(*func)(b,func2n),fc,p,q,rr,s,tol1,xm;
//        void nrerror(char error_text[]);

	if (fb*fa > 0.0) nrerror("Root must be bracketed in ZBRENT");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (fabs(fc) < fabs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*fabs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (fabs(xm) <= tol1 || fb == 0.0) return b;
		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				rr=fb/fc;
				p=s*(2.0*xm*q*(q-rr)-(b-a)*(rr-1.0));
				q=(q-1.0)*(rr-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=fabs(p);
			min1=3.0*xm*q-fabs(tol1*q);
			min2=fabs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (fabs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
		fb=(*func)(b,func2n);
	}
	nrerror("Maximum number of iterations exceeded in ZBRENT");
        return -1000.0;  // indicates error
}


#undef ITMAX
#undef EPS

void cajunctn(int t[],int ix,float r[5][3],float m[3][3],float det,float ang,
     float u1[3],float v1[3],float p1[3],
     float u2[3],float v2[3],float p2[3],float w1[3],float w2[3])
/*
NEW VERSION: arbitrary variable assigned = x - coef = b2, 
     needed if determinant = 0
r[0] = u1, r[1] = v1, r[2] = -u2  r[3] = -v2
 Keep old u1du2 etc. definitions but must change some signs 
   becuase a wnat r1 dot r2 etc. 
 
 routine to find the Ca position connnects two adjacent
  peptide planes.
Given vectors  u1 and v1 in plane 1 and u2 , v2 in plane 2
  v1 and v2 end at point p1,p2.
  find  w1 and w2  such that p1+w1 = Ca = p2 + w2
Using fact that w1 must be in u1,v1 plane: w1 = a u1 + b v1
   and w2 in u2,v2 plane and angle (w1 dot w2) must = ang;
 Can uniquely find Ca.
Procedure:  4 varialbles,  a1,b1,a2,b2,
   Have 3 linear equations and 1 non-linear
  Solve the 3 linear equaatons for a1,b1,a2 as function of b2
  then Insert these values into 4th eq. and find root (x = b2)
    using zbrent. (should be bounded by 0< b2 < 2)
*/
{
   int i;
   float cosang,b2;
   float  u1du1,u1dv1,v1dv1;  // u1 dot u1 etc.
   float  u2du2,u2dv2,v2dv2;  // u1 dot u1 etc.
   float u1du2,u1dv2,v1du2,v1dv2,dp[3],tol;
   float cc[4];  // the new general coefficients

   tol=1.0e-4;
   u2du2 = dot(u2,u2);
   u2dv2 = dot(u2,v2);
   v2dv2 = dot(v2,v2);
   u1du1 = dot(u1,u1);
   u1dv1 = dot(u1,v1);
   v1dv1 = dot(v1,v1);
   u1du2 = dot(u1,u2);
   u1dv2 = dot(u1,v2);
   v1du2 = dot(v1,u2);
   v1dv2 = dot(v1,v2);
   for(i = 0;i<=2;i++)
     dp[i] = p2[i] - p1[i];
   cosang = cos(ang);
#if 0
   for(i= 0;i<=100;i++){
      float xi = (i-50.0)/25.0;
      fprintf(fpout,"b2func( %f) = %f\n",xi,b2funcn(xi,func2n) );
   }
#endif
   // Find bounding points - Know answer must be  >0 and < 2, start at 0, step 0.1

   float x0 = 0.0;
   float xstep = 0.05;
   float y0 = b2funcn(x0,func2n);
   float y1 = b2funcn(x0+xstep,func2n);
//fprintf(fpout,"y0 = %f  y1 = %f\n",y0,y1);
   if(fabs(y0)<tol){
      b2 = y0;
      goto skipzbrent;
   }
   if(fabs(y1)<tol){
      b2 = y1;
      goto skipzbrent;
   }
   int isign,iplus;
   float xlow,xhigh;
   if(y0<y1){
       if(y0<0){ // increase  x1 in steps till sign changes
         isign =1;
         iplus= -1;
       }
       else{ // decrease  x1 in steps till sign changes
         isign = -1;
         iplus=1;
       }
   }
   else{ // y0>y1
       if(y0<0){ // increase  x1 in steps till sign changes
         isign = -1;
         iplus= -1;
       }
       else{ // decrase  x1 in steps till sign changes
         isign = 1;
         iplus=1;
       }
   }
fprintf(fpout,"isign = %d  iplus = %d\n",isign,iplus);
   float x1,x,y;
   x1 =x0;
   y1 = y0;
   do{// when exit this loop, have changed sign or have reached minmum value and now exact solution
        // first check for change of sign, then for reaching minimjm
     y = y1;  // store previous value
     x = x1;
     x1 = x1+isign*xstep;
     y1 = b2funcn(x1,func2n);
     if(x1>2.0){
       fprintf(fpout,"ERROR: x1 = %f CANNOT FIND 0 FOR zbrent in cajunct\n",x1);
       exit(0);
     }
     if(x1< -2.0){
       fprintf(fpout,"ERROR:x1 = %f CANNOT FIND 0 FOR zbrent in cajunct\n",x1);
       exit(0);
     }
     if(iplus*y1<0){ // then sign has changed  bounded by x,x1
        goto dozbrent;
     }
     if(fabs(y1)>fabs(y)){ // then has reached minimum - no exact solutin
fprintf(fpout," skipping zbrent becuase no exact solution\n");
       b2 = x; // Previous value - which must be minimum
       goto skipzbrent;
     } 
   }while(1);  // exit loop only with goto statments
   dozbrent:;
   if(isign>0){
       xlow = x;
       xhigh = x1;
   }
   else{
       xlow = x1; 
       xhigh = x;
   }
   b2=zbrentn(b2funcn,xlow,xhigh,func2n,tol);
   skipzbrent:;  // skip to here if use approximate method initial quess is corect answer
//WARNING: FOLLOWING PRINT STATEMENT CALLs b2func - sets values of cc
fprintf(fpout," At end of zbrent: cc[0]= %f cc[1] = %f cc[2] = %f cc[3] = %f  value = %f\n",
       cc[0],cc[1],cc[2],cc[3],b2funcn(b2,func2n) );
   for(i=0;i<=2;i++){
      w1[i] = cc[0]*r[0][i]+cc[1]*r[1][i];
      w2[i] =  -cc[2]*r[2][i]-cc[3]*r[3][i];
   }
#if 0 // check answer
float chpos[3];
   for(i=0;i<=2;i++)
      chpos[i] = p1[i]+w1[i];
   fprintf(fpout,"Ca theory (w1) pos =  %f %f %f  \n",chpos[0],chpos[1],chpos[2]);
   float mp[3];
   for(i=0;i<=2;i++)
     mp[i] =(p1[i]+p2[i])/2.0;
   fprintf(fpout,"Mid point:  %f %f %f\n",mp[0],mp[1],mp[2]);
   fprintf(fpout,"First point:  %f %f %f\n",p1[0],p1[1],p1[2]);
//   for(i=0;i<=2;i++)
//      chpos[i] = p2[i]+w2[i];
//   fprintf(fpout,"Ca theory (w2) pos =  %f %f %f  \n",chpos[0],chpos[1],chpos[2]);
#endif
}

 
 void findphipsierr(struct griddata& den1,struct geometry& geo)
/* routine to estimate phi/psi from bone pts.
Procdure:
  1. Find bone pts in residue plane:
      Start with pt1,pt,pt  - make co-ordinate system as follows              
   Define av = unit vector from first to second pt.
          bv    ..              second to third
  Set up coordinate system with x = av; z = x cross bv (make unit lenght)
         and y = z cross a
  Then go to next pt - 
         cv = unit vector for third to forth.
         theta = ang between z and cv = acos( dot(cv,z) );
  If theta close to PI/2.0, then four pts from plane.
 2. Then keep adding pts until - vector not in plane (i.e. theta not = PI/2)
 3. Then the next 4 points should form another plane  - check. 
Now if given three points in each plane can:
  1) estimate phi, psi by assuming that vector just before Ca is
    algned along N - Ca and vector after Ca is aligned along Ca - C
    This work OK if vectors are close to ideal, however can get large error
      if use big step  ( see routine findchivec below)
  2) Can find exact position of Ca between 2 vectors in each plane before
       and after Ca.  Based on idea that if p1 and p2 are points just
       before and after Ca, and w1 and w2 are vectors from p1 to Ca and
       p2 to Ca, just have
       p1 + w1 = p2 +w2 = Ca  -   This provides 3 equations.
       also have fourth equation   w1 dot w2 = N - Ca - c angle
    Can solve for 2 unknowns:  w1 = a1 u1 + b1 v1; w2 = a2 u1 + b2 v2,
      where u,v are 2 vectors in each plane and w must also be in this
      same plane, so onlh  2 unknowns for each plane.
      (see my routine  cajunct);
  Define pos[i] 1=ca,2=c,3=o;4=n,5=cb; 6=cb2
for i: phi is angle formed by Ci-Ni+1-Cai+1-Ci+1
       psi is angle formed by Ni+1-Cai+1-Ci+1-Ni+2
{
   findchi(pos[2][i],pos[4][i+1],pos[1][i+1],pos[2][i+1],phi);
   findchi(pos[4][i+1],pos[1][i+1],pos[2][i+1],pos[4][i+2],psi);
}
   chi is defined by 3 vectors:
  a) for phi:  vector: u= Ci to Ni+1,v= Ni+1 to Cai+1 and w= Cai+1 to Ci
    where Cai+1 to Ci is vector in plane after plane change.
    (I think any vector in this plane will give same phi)
    But, must have correct vector v= Ni+1 to Cai+1 --
     This is where error enters - since, depending length of step and
       its position, the last vector before plane change will not
       be exactly correct.
  b) for psi: 
        u = Ni+1 to Cai+1, v = Cai+1 to Ci+1, and w = Ci+1 to Ni+2
    (I think any vector in first plane can be used for u, and any
       vector second plane for w, but, again, must have correct v
      this is were error produced. 
 For this routine:
   Input is six points  p0 to p6.
    p0,p1,p2 are points in first plane, closest to theta change
    p3,p4,p5 are points in second plane closest to theta change

Procedure to find range of error in phi/psi
 First as check:
  Starting at Ca2 ( position of change of plane),
   go back 3 points and forward 2 points
  these 5 points should give  perfect phi/psi.
  The move back varying distance from Ca as start- and then go back 3 points,
    then forward 3 points, and use last 2 points to determine phi/psi.   
*/
{ 
   int i,j,k,ii,ii0,ivec[3],p,st[3],atp,resp,atn,resn,resf,af2;
   int p1,ix,t[3];
   float fvec0[3],fvec[3],len[10];
   float spos[10][3];
   float phi,psi,phirev,psirev,uv[6][3],cav[3],errdist,dist,ang,dist1,cv[3];
   float u1[3],u2[3],v1[3],v2[3],p0[3],p2[3],w1[3],w2[3],phi0,psi0;
   float m[3][3],det,r[5][3];

   fprintf(fpout," Finding error in phi/psi estimate\n");
   ivec[0]=ivec[1]=ivec[2]=0;
   imcv(den1,ivec,fvec0);// convert int grid ivec to float world co-ord fvec
   for(i=0;i<=1;i++) // find lengths as function of p
     for(j=0;j<=1;j++)
       for(k=0;k<=1;k++){
         p=i +2*j +4*k -1;
         if(p>=0){
            ivec[0]=i;ivec[1]=j;ivec[2]=k;
            imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
            len[p] = distv(fvec0,fvec);
         }
       }
    st[0]=1; // Ca positon  Forward direction : 0 to 1 to 2  (Ca to C to N of next residue)
    st[1]=2; // C 
    st[2] = 4; // N
// Set up test start positons at Ca, C and N and intermediate positons
    int resp0=3; // The residue number of the prev atom position
    for(i=0;i<=2;i++){ 
        cav[i] = -geo.pos[resp0][1 ][i] +geo.pos[resp0][4 ][i];// make unit vector from N to Ca
        cv[i] = geo.pos[resp0][2 ][i]-geo.pos[resp0][1 ][i];
    }
    dist = sqrt( dot(cav,cav) );
    dist1 = sqrt( dot(cv,cv) );
    ang = acos( dot(cav,cv)/(dist1*dist) );
fprintf(fpout,"ang = %f tt[4] = %f\n",ang,geo.tt[4]); 
    for(i=0;i<=2;i++)
         cav[i] = cav[i]/dist;
    errdist = 0.0; // distance in angstroms of start position before Ca
errdist = 0.5; // distance in angstroms of start position before Ca
    for(i=0;i<=2;i++) // starting point
      spos[3][i] = geo.pos[resp0][1 ][i]+ errdist*cav[i]; // Positiion near Ca -corresponds to point #3
    atp= 2; // The prev atom number  (N)
    resp = resp0;
    p1 = 6;
    p1 = 0;
    p1 = 6;
    ii0=ii=3;
    findprevpos(geo,resp,atp,&resf,&af2,st,len[p1],spos[ii],spos[ii-1]);
fprintf(fpout,"after first previous  res = %d  atom = %d\n",resf,af2);
// check proceudre:
    atp = af2;
    resp = resf;
    findprevpos(geo,resp,atp,&resf,&af2,st,len[p1],spos[ii-1],spos[ii-2]);
    atp = af2;
    resp = resf;
fprintf(fpout,"after second previous  res = %d  atom = %d\n",resf,af2);
    findprevpos(geo,resp,atp,&resf,&af2,st,len[p1],spos[ii-2],spos[ii-3]);
#if 0 // case for errdist = 0
    atn = 1;
#endif
#if 1 // case for errdist not = 0
    atn = 0;
#endif
    ii = ii0;    
    resn =resp0;
    findnextpos(geo,resn,atn,&resf,&af2,st,len[p1],spos[ii],spos[ii+1]);
fprintf(fpout,"Ca to first next = %f\n",distv(geo.pos[resp0][1 ],spos[ii+1]) );
    atn = af2;
    resn = resf;
    findnextpos(geo,resn,atn,&resf,&af2,st,len[p1],spos[ii+1],spos[ii+2]);
fprintf(fpout,"after second next  res = %d  atom = %d\n",resf,af2);
    atn = af2;
    resn = resf;
    findnextpos(geo,resn,atn,&resf,&af2,st,len[p1],spos[ii+2],spos[ii+3]);
fprintf(fpout,"after third next  res = %d  atom = %d\n",resf,af2);
    atn = af2;
    resn = resf;
    findnextpos(geo,resn,atn,&resf,&af2,st,len[p1],spos[ii+3],spos[ii+4]);
fprintf(fpout,"after fourth next  res = %d  atom = %d\n",resf,af2);
   for(i=0;i<=6;i++)
     for(k=0;k<=2;k++)
       uv[i][k] = spos[i+1][k]-spos[i][k];

// check if vectors in plane
float perpv[3][3],cosang;
   crossprod(cv,uv[5],perpv[0]);
   dist = sqrt(dot(perpv[0],perpv[0]) );
   for(i=0;i<=2;i++)// unitize perpv
      perpv[0][i] = perpv[0][i]/dist;
   dist= sqrt( dot(uv[6],uv[6]) );
   cosang = dot(uv[6],perpv[0])/dist;
fprintf(fpout,"cosang = (deg) = %f\n",cosang*180/PI);


//DEFINE PARAMETERS USED TO INTERPOLATE TO FIND Ca
   ang= geo.tt[4]; // C - Ca - N angle
   for(i=0;i<3;i++){
     u1[i] = uv[1][i];
     v1[i] = uv[2][i];
     u2[i] = -uv[5][i];
     v2[i] = -uv[4][i];
     p0[i] = spos[3][i];
     p2[i] = spos[4][i];
   }
#if 0 // test reverse direction
   for(i=0;i<3;i++){
     u1[i] = -uv[5][i];
     v1[i] = -uv[4][i];
     u2[i] = uv[1][i];
     v2[i] = uv[2][i];
     p0[i] = spos[4][i];
     p2[i] = spos[3][i];
   }
#endif
#if 0 // old version - sets b2 = x
    for(j=0;j<=2;j++){ // define each column  m[i][j]
       m[j][0] = u1[j];
       m[j][1] = v1[j];
       m[j][2] = -u2[j];
    }     
    det = determ(m);
fprintf(fpout," findphi det = %f\n",det);
    if(fabs(det)< 0.1)
       fprintf(fpout," ERROR: det in findphipsierr = %f < %f\n",det,0.1);
   cajunct(m,det,ang,u1,v1,p0,u2,v2,p2,w1,w2);
#endif

#if 1 // New Version - choose x arbitrarily
   for(i=0;i<=2;i++){  // define r vectors 
      r[0][i] = u1[i];
      r[1][i]= v1[i];
      r[2][i] = -u2[i];
      r[3][i] = -v2[i];
   }
for(ix = 0;ix<=3;ix++){
//for(ix = 0;ix<=0;ix++){
   k= -1;
   for(i=0;i<=2;i++){// set 3 dependent variable = t[0,1,2]
     if(i==ix)
       k = k+2;
     else
       k = k+1;
     t[i] = k;
   }
   for(j=0;j<=2;j++){ // define each column  m[i][j]
     for(i= 0;i<=2;i++)
       m[j][i] = r[ t[i] ][j];
   }
    det = determ(m);
fprintf(fpout," findphi:  ix = %d  det = %f\n",ix,det);
   cajunctn(t,ix,r,m,det,ang,u1,v1,p0,u2,v2,p2,w1,w2);
} // ix
#endif


   phi0= geo.phihel*180.0/PI;  // theoretical hellix values
   psi0= geo.psihel*180.0/PI;
fprintf(fpout," Actual Ca pos = %f %f %f\n",geo.pos[resp0][1 ][0],geo.pos[resp0][1 ][1],geo.pos[resp0][1 ][2]);
#if 0 // used for errdist = 0
   findchivec(uv[1],uv[2],uv[3],&phi);
   findchivec(uv[2],uv[3],uv[4],&psi);
#endif
#if 0 // used for errdist not 0  (change in plane across vector 3)
   findchivec(uv[1],uv[2],uv[4],&phi);
   findchivec(uv[2],uv[4],uv[5],&psi);
 #endif
#if 1 // Use interpolated CA (change in plane across vector 3)
   for(i=0;i<=2;i++)
     w2[i] = -w2[i]; // vector from Ca top pos 4
//   findchivec(uv[2],w1,w2,&phi);
   findchivec(uv[1],w1,w2,&phi);  // use this if uv[2] = w1
   findchivec(w1,w2,uv[4],&psi);
#endif
   fprintf(fpout,"\n\n USING CA:  step length = %f\n",len[p1]);
   fprintf(fpout,"step length = %5.3f erdist = %5.3f:  phi = %5.3f  psi = %5.3f (should be %f %f)\n",
      len[p1],errdist,phi*180.0/PI,psi*180.0/PI,phi0,psi0);
   for(i=0;i<=5;i++) // find phi psi for the reverse direction
     for(k=0;k<=2;k++)
       uv[i][k] = -uv[i][k]; // change sign of all vectors
#if 1 // used for errdist not 0  (change in plane across vector 3)
   findchivec(uv[5],uv[4],uv[2],&phirev);
   findchivec(uv[4],uv[2],uv[1],&psirev);
#endif
   fprintf(fpout,"\n\n MAIN RESULTS:  step length = %f\n",len[p1]);
   fprintf(fpout,"step length = %5.3f erdist = %5.3f:  phi = %5.3f  psi = %5.3f phirev= %5.3f psirev = %5.3f (should be %f %f)\n",
      len[p1],errdist,phi*180.0/PI,psi*180.0/PI,phirev*180.0/PI,psirev*180.0/PI,phi0,psi0);
#if 1 // Use interpolated CA (change in plane across vector 3)
   for(i=0;i<=2;i++)
     w2[i] = -w2[i]; // vector from Ca top pos 4
   findchivec(uv[2],w1,w2,&phi);
   findchivec(w1,w2,uv[4],&psi);
#endif
   fprintf(fpout,"\n\n USING CA:  step length = %f\n",len[p1]);
   fprintf(fpout,"step length = %5.3f erdist = %5.3f:  phi = %5.3f  psi = %5.3f phirev= %5.3f psirev = %5.3f (should be %f %f)\n",
      len[p1],errdist,phi*180.0/PI,psi*180.0/PI,phirev*180.0/PI,psirev*180.0/PI,phi0,psi0);
}

void nclassifyphipsi(struct geometry& geo,float phi,float psi,
   int *type)
/*
MODIFY  from classify - 11aug99, if nclass = 0 set type = 0
  Used in writepdbtofit
 */
{
   int i;

   for(i=1;i<=geo.ntypes;i++){
     if( (phi> geo.sheet[i][0])&&(phi<geo.sheet[i][1])&&
               (psi> geo.sheet[i][2])&&(psi<geo.sheet[i][3]) ){
        *type = i;
        goto quitrout;
     }
   }
   *type = 0; // if reach here - not in any class
   quitrout:;
}

void classifyphipsi(struct geometry& geo,float phi,float psi,
   int *nclass,int type[],float radius[])
/*
 Classify phi/psi angle into constraint type:
Type 1 = helix, 2 = sheet, 3  = left helix 4= epsilon 5 = delta  6= gammar 7 = gammal 
Returns: nclass = number of constraints for which angles fall within min,max 
                     range
         type[1,2,..nclass] = types for each of the above classes
         radius[1,2,3,..geo.ntypes] = the radial distance from center of each type
    IF nclass = 0 , type[1,2] = type of closest and next closest region
    If nclass = 1;  type[1] = type for nclass, type[2] = next closest radiee
*/
{
   int i,ttype[4];
   float minrad1,minrad2;

   *nclass = 0;
   minrad1 = minrad2 = 10000.0;
   for(i=1;i<=geo.ntypes;i++){
     radius[i] = sqrt( sqr(phi -geo.phi[i])+sqr(psi-geo.psi[i]) );
     if( (phi> geo.sheet[i][0])&&(phi<geo.sheet[i][1])&&
               (psi> geo.sheet[i][2])&&(psi<geo.sheet[i][3]) ){
        *nclass = *nclass+1;
        type[*nclass] = i;
     }
   }
   if(*nclass <=1){ // find two closest radius
     for(i=1;i<=geo.ntypes;i++){
        if(radius[i]< minrad1){ // minrad1 = closest radius
           if(i>1){ // set second closest = previus first closest
             ttype[2] = ttype[1];
             minrad2 = minrad1;
           }
           ttype[1]= i;
           minrad1 = radius[i];
        }
        else if(radius[i]< minrad2){// minrad2 = next closest radius
           ttype[2] = i;
           minrad2 = radius[i];
        }
      }
      if(*nclass==1){ // then set type[2] = min radius not in min/max region
         if(ttype[1]!= type[1])
             type[2] = ttype[1];
         else
             type[2] = ttype[2];
      }
      if(*nclass==0){
          type[1] = ttype[1];
          type[2] = ttype[2];
      }
    }
}
           
     
