// gl version of thinpdbtolev    use minimum of global variable
#if 0
Main definition of mainch and sidech:
ii = segment number    for(ii=0;ii<=nseg;ii++)
i = res number          for(i=1;i<=totres[ii];i++)
    for(j=0;j<=2;j++)
        mainch[ii][i][0][j]=res[ii][i].pos[N][j];
    for(j=0;j<=2;j++)
        mainch[ii][i][1][j]=res[ii][i].pos[CA][j];
    for(j=0;j<=2;j++)
       mainch[ii][i][2][j]=res[ii][i].pos[C][j];


   resdat[ii+1][i][0]=nbrn+1;
      or nbrn = resdat[ii+1][i][0]-1
   for(j=1;j<=nbrn;j++) { /*jth branch*/
          npts=side[iresname][j][0];
          resdat[ii+1][i][j]=npts;
          for(k=1;k<=npts;k++) {/* kth point*/
             atomname=side[iresname][j][k];             
             for(jj=0;jj<=2;jj++){ /* x,y,z */
                sidech[ii][i][j][k-1][jj]=res[ii][i].pos[atomname][jj];
    j=nbrn+1
     for(jj=0;jj<=2;jj++) 
        sidech[ii][i][j][0][jj]=res[ii][i].pos[C][jj];
     for(jj=0;jj<=2;jj++) 
        sidech[ii][i][j][1][jj]=res[ii][i].pos[O][jj];
#endif

/*
NEW:  21jun00 - USE TER as indication of new segment - just
   ignore segment id.

NOTE: THERE IS SOME KIND OF ERROR IN LINE BELOW (IN newpdbtolev):
                    subst[nsubst].pos[0][0]=atomnum;
   1) sets integer to float.
   2) It is necessary.  It seems to work OK. Have not changed it. 

New version  - does not make levfile, but does everything in getmodel
   New to comment out end of getmodel

New version used as run time program by thin

New version - reads either xplor or brookhaven pdb files
New Version - prints out actual residue number at end of first line
   for residue.  Used in maid to print actual number
 NEW version - 
1)finds gaps in chain 
   All atoms whose separation
     is less than "bondlengh" (2.0 A) are assumed to bonded
2) Remove the "x" and "s" input, use only xplor format
3) Comment out printer "WARNING" if numbers are skipped
Note: OT2 (terminal carboxyl O) is not drawn

NOTE: Corrected error in version now installed: should be  <=2 for min, max[i]
  Procedure for adding additional substrates:
    1) Add the new 3 letter substrate name to residuenum list, after amino 
       and before SOL.
    2) Add the new atom names to atomnum list, in any posiiton,but try to 
       keep alphabetical.
    3) Add the substrate structure to the data array as follows:
       a)name; b) # of branchs (treat entire structure as branches)
       c)# of atoms in first branch; d)atom list of first branch in order;
       e)# of atoms in second branch; etc.
       Make the first branch the dominant chain of substrate, the center of
         this branch will be used for positioning substrate name on screen
     4) Add the necessary "else" statements to assignatomnum and assignresnum
        routines.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pdbtolev.h"
#include "drawwindow.h"
#include "math.h"
#include "peptide.h"
#include "fittopdb.h"

 extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version

int side[AA][NBRN][NPTS];  // defined globally because needed in assignseg.c++ 
int brdraw[AA][NBRN]; // defined globally because MAYBE?? needed in assignseg.c++

void bond(struct substrate* subst,int side[AA][NBRN][NPTS],int snum,int subname,int natot,float bsqr,int i,int i1,int *bi); 
void myconnect(struct substrate* subst,int side[AA][NBRN][NPTS],int snum,int subname); /*snum=substrate number and subname=101*/
#if 0 // listd in pdbtolev.h
void checkatomname(int side[AA][NBRN][NPTS],char nres[],int resname,int *atomnum);
void assignresnum(char nres[],char nm[],int *resn);
void assignatomnum(char nres[],char resname[],char na[],int *atomn);
#endif

const char*  nameres[30] = {"ZERORES","ALA","ARG","ASN","ASP","CYS","GLU","GLN","GLY","HIS","ILE",
          "LEU","LYS","MET","PHE","PRO","SER","THR","TRP","TYR","VAL","AA3","AA4"};

const char* nameatom[100] = { "ZER0","CA","C","O","N","CB","OG","CD","CD1","CD2","CE","CE1","CE2","CE3",
           "CG","CG1","CH1","CH2","CG2","CZ","CZ1","CZ2","CZ3","C1","C2","C3",
           "C4","C5","C6","ND1","ND2","NE","NE1","NE2","NH1","NH2","NZ","OD","OD1",
           "OD2","OE1","OE2","OG1","OG2","OH",
           "O11","O12","O3","O41","O42","O61","O62","SD","SG","OHH" };



void conv_lev(struct residue **res,struct substrate *subst,
        float ****subch,char ***rname,
        char** sname,int **subdata,char ***resnumch,
         float ****mainch,float *****sidech,
         int **gaps,
         char** segname,
        int ***resdat,float **maincolor,float **sidecolor,
        int side[AA][NBRN][NPTS],
        int nseg,int totres[],int totsubst,
        char** namelist,char** subnamelist); /* convert to lev format*/
float tsqr(float x); /* given different name so does not conflict with previously defined sqr*/

struct residue 
{
  int num; /* enum residue number, see below*/
  char actnum[5]; /* the real number of the residue as a character array*/
  float pos[NATOMS][3]; /*x,y,z position of atom using enum atomnum*/
}; 

struct substrate
{
  int num;
  float pos[NATOMS][3];
}; 
 
//int data[AA+NSUB][50]={   Original version  - do not add substrate to dat ??
int data[AA+4][50]={
     {ALA,1,2,CA,CB},
     {ARG,2,7,CA,CB,CG,CD,NE,CZ,NH1,2,CZ,NH2},
     {ASN,2,4,CA,CB,CG,OD1,2,CG,ND2},
     {ASP,2,4,CA,CB,CG,OD1,2,CG,OD2},
     {CYS,1,3,CA,CB,SG},
     {GLU,2,5,CA,CB,CG,CD,OE1,2,CD,OE2},
     {GLN,2,5,CA,CB,CG,CD,OE1,2,CD,NE2},
     {GLY,0},
     {HIS,1,8,CA,CB,CG,ND1,CE1,NE2,CD2,CG},
     {ILE,2,4,CA,CB,CG1,CD1,2,CB,CG2},
     {LEU,2,4,CA,CB,CG,CD1,2,CG,CD2},
     {LYS,1,6,CA,CB,CG,CD,CE,NZ},
     {MET,1,5,CA,CB,CG,SD,CE},
     {PHE,1,9,CA,CB,CG,CD1,CE1,CZ,CE2,CD2,CG},
     {PRO,1,5,CA,CB,CG,CD,N},
     {SER,1,3,CA,CB,OG},
     {THR,2,3,CA,CB,OG1,2,CB,CG2},
     {TRP,2,12,CA,CB,CG,CD1,NE1,CE2,CZ2,CH2,CZ3,CE3,CD2,CG,2,CD2,CE2},
     {TYR,2,7,CA,CB,CG,CD1,CE1,CZ,OH,4,CG,CD2,CE2,CZ},
     {VAL,2,3,CA,CB,CG1,2,CB,CG2},
     {CIT,6,7,O11,C1,C2,C3,C5,C6,O61,2,C1,O12,3,C3,C4,O41,2,C4,O42,2,C3,O3,2,C6,O62}
}; 


float bondlength;
int cmat[NATOMS][5];

void resmat2(struct residue** &mat,int n1,int n2)// allocates array 2 dim array (mat)
         // of type int - n1,n2,n3,n4 size of each dimension
{
   int i;
   mat =new residue*[n1];
   assert(mat  !=0);
   for ( i = 0; i < n1; ++i){
      mat [i] = new residue[n2];
      assert(mat [i] != 0);
   }
}


void rereadpdb(struct pdbfile& pdb1)
/* 
   Just runs pdbtolev again - assumes that
    no change in residue number, etc. from when allocatelev was run
*/
{
     int i,nseg,nres[1000],nsubst,nwater,maxres; 
     findpdbdimensions(pdb1.name,&nseg,nres,&nsubst,&nwater);
     for(i=0;i<nseg;i++)
        if(nres[i]>maxres)
          maxres= nres[i]+2;
      pdbtolev(nseg,nwater,maxres,nsubst,pdb1.subch,pdb1.rname,pdb1.subdata,pdb1.resnumch,
           pdb1.mainch,pdb1.sidech,pdb1.gaps,pdb1.resdat,
           pdb1.maincolor,pdb1.sidecolor,pdb1.name,pdb1.xmin,pdb1.xmax);
}


void pdbtolev(int countseg,int countwater,int countres,int countsubst,
    float ****subch,char ***rname,
    int **subdata,char ***resnumch,float ****mainch,
    float *****sidech,int **gaps,int ***resdat,
    float **maincolor,float **sidecolor,
       char pdbfile[],float xmin[3],float xmax[3])
/* NEW VERsion - 21jun00 - define segments by presence of "TER" 
    statement - complete ignore segment id
  Only change segment number if have read at least 1 ATOM before ter
    Thus, will not change for substrates, etc.
*/
{
   int readatom,i,nseg,nres,nsubst,atomnum;
   int *totres = new int[countseg];
   int nwater,rnum,rtype,totsubst;
//   int brdraw[AA][NBRN];// Now defined globally
   char** sname;
   char** segname;
   FILE *fp;
   struct residue** res;
   struct substrate* subst;
   float x[3],sigma,bfact;
   float** xwater;
   char lresnumch[NCHAR],oldresnumch[NCHAR]; 
   char dataline[200],atomch[30],atomnumch[NCHAR],atomname[NCHAR],
       resname[NCHAR],seg[NCHAR],oldseg[NCHAR];
   char** namelist;
   char** subnamelist;

   cmat2(segname,countseg+1,20); //Allocate memory
printf(" countres = %d, NCHAR = %d\n",countres,NCHAR);
   cmat2(namelist,VAL+1,NCHAR);
   fmat2(xwater,countwater+1,3);
   resmat2(res,countseg,countres+1);
   if(countsubst>0){ //NOTE THIS ALLOCATION IS WRONG, SEE namelist
     subst = new substrate[countsubst+1];
     cmat2(subnamelist,countsubst+1,NCHAR);   
     cmat2(sname,countsubst+1,4);
   }

 bondlength=2.0; /*distance in Angstroms for which 2 atoms are considered bonded*/
   if( (fp=fopen(pdbfile,"r") ) ==NULL) {
      fprintf(stderr,"error opening file %s\n",pdbfile);
      exit(0);
   }
 /* convert from data[RES][list#] to side[RESNAME][..][..] where
     side[res name][0][0]=# of branches in  residue name;
     side[name][i>0][0]=# of points in ith branch of res neame;
     side[name][i>0][j>0]= atom name of jth point in ith branch of res name
  */ 
/*  */ /* now done at beginning at thinchain.c */
  data_side(side,brdraw);
   /* read into res[nseg][1... total # of residues]*/
   strcpy(oldresnumch,"");
   nres=0;
   nsubst=0;
   for(i=0;i<=2;i++) {
       xmin[i]= 1000.0;
       xmax[i]= -1000.0;
   }
   nseg= -1; // first segment number (0,1,2..
//   itoa(nseg+1,segname[nseg]); // COMMENTED OUT ON 13nov00- seems to be big error
   nres=0;
   nwater=0;
   strcpy(oldseg,"");
   strcpy(seg,"");/* need this because if seg part of data line missing, will use whatever is intial value*/
   readatom = 0;
   while(  fgets(dataline,199,fp) != NULL) {
    if(strncmp("TER",dataline,3)==0){ // If line starts with TER, new segment
       readatom = 0;
    }
    else if(strncmp("ATOM",dataline,4)==0){ /*Read only if begins with ATOM*/
        if(readatom==0){
           readatom = 1;
           if(nseg != -1)
                   totres[nseg]=nres; /* nres from previous read*/
           nseg = nseg+1;
           nres = 0;
           itoa(nseg+1,segname[nseg]);
        }
        if(dataline[21]!=' ') /* then standard pdb format */
            sscanf(dataline,"%s%s%s%s%s%s%8f%8f%8f",atomch,atomnumch,atomname,
               resname,seg,lresnumch, &x[0],&x[1],&x[2]);
        else{ /* xplor format */
           sscanf(dataline,"%s%s%s%s%s%8f%8f%8f%6f%6f%s",atomch,atomnumch,atomname,resname,lresnumch,
                                     &x[0],&x[1],&x[2],&sigma,&bfact,seg);
       }
#if 0
       if( strcmp(seg,"") == 0) /* segname blank*/
           nseg=0;
       else { /*sement name not blank*/
            if( (seg[0]<'A')||(seg[0]>'Z') ) /* 65 is character A, if not A...Z then assume null */
                 nseg=0;
            else if( strcmp(seg,oldseg)!=0) {/*change in segment name,*/
               if (nres!=0)  /*not first line*/
                   totres[nseg]=nres; /* nres from previous read*/
               strcpy(oldseg,seg);
               nseg=nseg+1;
               strcpy(segname[nseg],seg);/* assign name of previous seg*/
               nres=0;        
             }
        }
#endif 
        /*discard H atoms and drop OT2,change OT1 to O (xplor format
             in pdb format, drop OT*/
       if( (strncmp(atomname,"H",1)!=0)&& (strcmp(atomname,"OT")!=0)&&
           (strcmp(atomname,"OT2") !=0)&&(strcmp(atomname,"OXT") !=0)&&
           (strncmp(atomname,"1H",2)!=0)&&(strncmp(atomname,"2H",2)!=0)&&
            (strncmp(atomname,"3H",2)!=0) ) {
         if((strcmp(resname,"SOL")==0)||(strcmp(resname,"WAT")==0)
                   ||(strcmp(resname,"HOH")==0) ){ /*if water,   just keep list of x,y,z pos*/
           
             nwater=nwater+1;
             for(i=0;i<=2;i++) 
               xwater[nwater][i]=x[i];
           }
         else {
           if(strcmp(atomname,"OT1") ==0)
                  strcpy(atomname,"O");
           if(strcmp(resname,"ILE")==0) /* in ILE, allow either CD or CD1*/
              if(strcmp(atomname,"CD")==0)
                 strcpy(atomname,"CD1");
           for(i=0;i<=2;i++) {
             if(x[i]<xmin[i]) xmin[i]=x[i];
             if(x[i]>xmax[i]) xmax[i]=x[i];
           } 
           if(strcmp(oldresnumch,lresnumch)!=0){ /*new residue*/
              assignresnum(lresnumch,resname,&rnum);/*if not in dictionary, 
                               assigned rnum=101*/
              if(rnum<=VAL){ /*amino acid, part of main chain*/
                rtype=0;/*rtype of 0 is amino acid*/
                nres=nres+1; 
                res[nseg][nres].num=rnum;
                strcpy(res[nseg][nres].actnum,lresnumch);
                /*name list is character names for each number*/
                for(i=0;i<=2;i++) /*read 3 characters of residue name*/
                   namelist[res[nseg][nres].num][i]=resname[i];
              }
              else { /*not amino acid,if not in dictionary rnum=101 */
                rtype=1; /* i.e. a "substrate"*/
                nsubst=nsubst+1;
                subst[nsubst].num=rnum;
                for(i=0;i<=2;i++) /*read 3 characters of residue name*/
                   subnamelist[nsubst][i]=resname[i];
              }
              strcpy(oldresnumch,lresnumch);
              atomnum=0; /*used only if not in dictionary,otherwise given new value*/
            }
            if(rnum<100){ /* in dictionary*/
               assignatomnum(lresnumch,resname,atomname,&atomnum);
               /*if amino acid, check atom assignment and make correction to ILE*/
               if(rtype==0)
                  checkatomname(side,lresnumch,res[nseg][nres].num,&atomnum);
               else
                 checkatomname(side,lresnumch,subst[nsubst].num,&atomnum);
               if(rtype==0){
                 for(i=0;i<=2;i++) 
                   res[nseg][nres].pos[atomnum][i]=x[i];
               }
               else { /*for a substrate*/
                 for(i=0;i<=2;i++) 
                    subst[nsubst].pos[atomnum][i]=x[i];
               }
             }
             else{ /*not in dictionary*/
                atomnum=atomnum+1; /*this was set to 0 at beginning of substrate*/
                for(i=0;i<=2;i++){ 
                    subst[nsubst].pos[atomnum][i]=x[i];
                }
                subst[nsubst].pos[0][0]=atomnum; /*total # of atoms*/
             }
         }/*if not water*/
       } /* if not H */                
      } /* if ATOM*/
   }
//   if(nseg==0) strcpy(segname[0],"A"); /* default name if only one segment*/
   totres[nseg]=nres;  /* last residue*/
   totsubst=nsubst;
   conv_lev(res,subst,subch,rname,sname,subdata,resnumch,mainch,
         sidech,gaps,segname,resdat,maincolor,sidecolor,side,
         nseg,totres,totsubst,namelist,subnamelist);
   fprintf(fpout,"minx = %f  maxx= %f\n",xmin[0],xmax[0]);
   fprintf(fpout,"miny = %f  maxy= %f\n",xmin[1],xmax[1]);
   fprintf(fpout,"minz = %f  maxz= %f\n",xmin[2],xmax[2]); 
   fclose(fp); 
}


void findpdbdimensions(char pdbfile[],int *nseg,int nres[],int *nsubst,
        int  *nwater)
/*
 NEW VERsion - 21jun00 - define segments by presence of "TER" 
    statement - complete ignore segment id
  Only change segment number if have read at least 1 ATOM before ter
    Thus, will not change for substrates, etc.
  segments are numbered  [0,1,2 ... nseg-1]
*/
{
   int rnum;
   FILE *fp;
   float x[3],fjunk;
   char lresnumch[NCHAR],oldresnumch[NCHAR]; 
   char dataline[200],atomch[30],atomnumch[NCHAR],atomname[NCHAR],resname[NCHAR],
       seg[NCHAR],oldseg[NCHAR];
   int readatom;

   if( (fp=fopen(pdbfile,"r") ) ==NULL) {
      fprintf(stderr,"error opening file %s\n",pdbfile);
      exit(0);
   }
   /* read into res[nseg][1... total # of residues]*/
   strcpy(oldresnumch,"");
   nres[0]=0;
   *nsubst=0;
   *nseg= -1;
   *nwater=0;
   strcpy(oldseg,"");
   strcpy(seg,"");/* need this because if seg part of data line missing, will use whatever is intial value*/
   readatom = 0;
   while(  fgets(dataline,199,fp) != NULL) {
    if(strncmp("TER",dataline,3)==0){ // If line starts with TER, new segment
       readatom = 0;
    }
    else if(strncmp("ATOM",dataline,4)==0){ /*Read only if begins with ATOM*/
        if(readatom==0){
           readatom = 1;
           *nseg = *nseg+1;
           nres[*nseg] = 0;
        }
        if(dataline[21]!=' ') /* then standard pdb format */
            sscanf(dataline,"%s%s%s%s%s%s%8f%8f%8f",atomch,atomnumch,atomname,
               resname,seg,lresnumch, &x[0],&x[1],&x[2]);
        else /* xplor format */
           sscanf(dataline,"%s%s%s%s%s%8f%8f%8f%6f%6f%s",atomch,atomnumch,atomname,resname,lresnumch,
                                     &x[0],&x[1],&x[2],&fjunk,&fjunk,seg);
#if 0
       if( strcmp(seg,"") == 0) /* segname blank*/
           *nseg=0;
       else { /*sement name not blank*/
            if( (seg[0]<'A')||(seg[0]>'Z') ){ /* 65 is character A, if not A...Z then assume null */
//                 fprintf(fpout,"NOTE:  segment name does not begin with A to Z -asumme only 1 segment\n");
                 *nseg = 0;
            }
            else if( strcmp(seg,oldseg)!=0) {/*change in segment name,*/
               strcpy(oldseg,seg);
               *nseg= *nseg+1;
               nres[*nseg]=0;        
             }
        }
#endif 
        /*discard H atoms and drop OT2,change OT1 to O (xplor format
             in pdb format, drop OT*/
       if( (strncmp(atomname,"H",1)!=0)&& (strcmp(atomname,"OT")!=0)&&
           (strcmp(atomname,"OT2") !=0)&&(strcmp(atomname,"OXT") !=0)&&
           (strncmp(atomname,"1H",2)!=0)&&(strncmp(atomname,"2H",2)!=0)&&
            (strncmp(atomname,"3H",2)!=0) ) {
         if((strcmp(resname,"SOL")==0)||(strcmp(resname,"WAT")==0)
                   ||(strcmp(resname,"HOH")==0) ){ /*if water,   just keep list of x,y,z pos*/
           
             *nwater= *nwater+1;
           }
         else {
           if(strcmp(oldresnumch,lresnumch)!=0){ /*new residue*/
              assignresnum(lresnumch,resname,&rnum);/*if not in dictionary, 
                               assigned rnum=101*/
              if(rnum<=VAL){ /*amino acid, part of main chain*/
                nres[*nseg]=nres[*nseg]+1; 
//fprintf(fpout,"nseg = %d  nres = %d\n",*nseg,nres[*nseg]);
              }
              else { /*not amino acid,if not in dictionary rnum=101 */
                *nsubst= *nsubst+1;
              }
              strcpy(oldresnumch,lresnumch);
            }
         }/*if not water*/
       } /* if not H */                
      } /* if ATOM*/
   }
   *nseg= *nseg+1;
   fclose(fp); 
}


void checkatomname(int side[AA][NBRN][NPTS],char nres[],int resname,int *atomnum)
//int *atomnum,resname;
//char nres[];
{
  int i,j,error,numpts,numbr,atomname;
  error=0;
  /*check main chain*/
  atomname= *atomnum;
  if (atomname==C) error=1;
  if (atomname==N) error=1;
  if (atomname==O) error=1;
  if (atomname==CA) error=1; 
  /*check side chain*/
  numbr=side[resname][0][0];
  for(i=1;i<=numbr;i++) {
     numpts=side[resname][i][0];
     for(j=1;j<=numpts;j++) 
        if(atomname == side[resname][i][j]) 
             error=1;  
  }
  if( error != 1) { /*change from pdb to xplor format for ILE*/
     if( (resname==ILE)&&(atomname==CD1)) *atomnum=CD; 
     else
      fprintf(fpout,"ERROR %d in res= %s, in atomname = %d\n",error,nres,atomname);   
   }
} 

void conv_lev(struct residue **res,struct substrate *subst,
        float ****subch,char ***rname,
        char** sname,int **subdata,char ***resnumch,
         float ****mainch,float *****sidech,
         int **gaps,
         char** segname,
        int ***resdat,float **maincolor,float **sidecolor,
        int side[AA][NBRN][NPTS], // now defined globally
        int nseg,int totres[],int totsubst,
        char** namelist,char** subnamelist) /* convert to lev format*/
{
  int ii, i,j,jj,k,npts,nbrn,iresname,subname,atomname,nmain;
  int kk;
  char rrname[4];
  float distsq,dist,xsmax[3],xsmin[3];
  int countgaps, ngaps[500]; /* maximum of 500 gaps per chain */
  /* in following arrays - order as: residue name,number of side chains, #of points 
        in first chain, atoms in first chain;  # of points in second chain, etc*/
  /*outputs line file type .lev in following format:
      first line: number of points on main chain; # of segments; 
      second: segment name 1,  number of residues (=nres) in segment 1;name of seg 2; ete
      third: through nres = x,y,z of N,CA,C - mainchain
      then:  branch chains: #points, x,y,z of points
  */
 distsq=9.0; /* square of max allowed distance between neighbboring residues to be connected*/
 nmain=3; /*number of points in main chain of each residue*/
 resdat[0][0][0]=nseg+1;
 for(i=0;i<=nseg;i++){ 
   resdat[i+1][0][1]=nmain;
   resdat[i+1][0][0]=totres[i];
 }
 for(i=nseg+1;i>=1;i--)
   strcpy(segname[i],segname[i-1]); /* renumber from 1 to nseg+1 */
#if 0 // contains screen stuff, now called in model.c++ after call to pdbtolev
 i=0;
   for(k=0;k<=2;k++) { /* first seg is yellow and blue */
     maincolor[i][k]=yellowv[k];
     sidecolor[i][k]=lightbluev[k];
   }
 for(i=1;i<nseg;i++){
   for(k=0;k<=2;k++) { /* all seg after first are in cyan and pink */
     maincolor[i][k]=cyanv[k];
     sidecolor[i][k]=pinkv[k];
   }
 }
#endif
 for(k=0;k<=2;k++) { /* set initial values of min max of x,y,z*/
      xsmax[k]= -1000.0;
      xsmin[k]= 1000.0;
 } 
 for(ii=0;ii<=nseg;ii++){ /* find gaps in chain */
    countgaps=1; /* always skip connection to first N */
    ngaps[1]=1; /* position of first gap */
    for(i=1;i<totres[ii];i++){ /* find the  C[i]-N[i+1]  distance */
       dist=0;
       for(k=0;k<3;k++)
         dist=dist+tsqr(res[ii][i].pos[C][k]-res[ii][i+1].pos[N][k]);
       if(dist>distsq){
         countgaps=countgaps+1;
         ngaps[countgaps]=i+1;
       }
    }
    gaps[ii][0]=countgaps;
    for(i=1;i<=countgaps;i++)
       gaps[ii][i]=ngaps[i];
 }/* ii */
 for(ii=0;ii<=nseg;ii++) {
   for(i=1;i<=totres[ii];i++) { 
#if 0 // Cannot have array = 100 - Also this is not used anymore
if(i==290)
  printf("\n");
    resdat[ii+1][i][100] = 0; // This indicates has not been modified
#endif
    iresname=res[ii][i].num;
    nbrn=side[iresname][0][0]; 
    resdat[ii+1][i][0]=nbrn+1;
    /*write out mainchain atoms*/
    for(j=0;j<=2;j++)
        mainch[ii][i][0][j]=res[ii][i].pos[N][j];
    for(j=0;j<=2;j++)
        mainch[ii][i][1][j]=res[ii][i].pos[CA][j];
    for(j=0;j<=2;j++)
       mainch[ii][i][2][j]=res[ii][i].pos[C][j];
    for(j=0;j<=2;j++) 
         rrname[j]=namelist[iresname][j];/*char name of residue*/
    rrname[3]=0; /* terminate with null */
    for(kk=0;kk<3;kk++){ /* 3 is the number of mainpts */
       for(j=0;j<=2;j++){
         if(mainch[ii][i][kk][j]>xsmax[j]) xsmax[j]=mainch[ii][i][kk][j];
         if(mainch[ii][i][kk][j]<xsmin[j]) xsmin[j]=mainch[ii][i][kk][j];
       }
    }/*kk = mainpts number*/
    strcpy(rname[ii][i],rrname);
    strcpy(resnumch[ii][i],res[ii][i].actnum);
  /*write out sidechains*/
    for(j=1;j<=nbrn;j++) { /*jth branch*/
          npts=side[iresname][j][0];
          resdat[ii+1][i][j]=npts;
          for(k=1;k<=npts;k++) {/* kth point*/
             atomname=side[iresname][j][k];             
             for(jj=0;jj<=2;jj++){ /* x,y,z */
                sidech[ii][i][j][k-1][jj]=res[ii][i].pos[atomname][jj];
                if(sidech[ii][i][j][k-1][jj]>xsmax[jj]) xsmax[jj]=sidech[ii][i][j][k-1][jj];
                if(sidech[ii][i][j][k-1][jj]<xsmin[jj]) xsmin[jj]=sidech[ii][i][j][k-1][jj];
            }/* jj */
          }/* k */
     }/* j */
     /* add C-O branch*/
     resdat[ii+1][i][j]=2;
     for(jj=0;jj<=2;jj++) 
        sidech[ii][i][j][0][jj]=res[ii][i].pos[C][jj];
     for(jj=0;jj<=2;jj++) 
        sidech[ii][i][j][1][jj]=res[ii][i].pos[O][jj];
   }/*for i*/
  } /*for ii*/
  subdata[0][0]=totsubst;
  for(i=1;i<=totsubst;i++){ /*write out substrate pos*/
     for(j=0;j<=2;j++) 
         rrname[j]=subnamelist[i][j];/*char name of residue*/
     subname=subst[i].num;
#if 1
     if(subname>100) /*not in dictionary*/
        myconnect(subst,side,i,subname); /*define connectivity based on distance between atoms*/
#endif
     nbrn=side[subname][0][0];/*treat substrate as series of sidechains*/
     strcpy(sname[i],rrname);
     subdata[i][0]=nbrn;
     for(j=1;j<=nbrn;j++) {
        npts=side[subname][j][0];
        subdata[i][j]=npts;
        for(k=1;k<=npts;k++) {
           atomname=side[subname][j][k]; /*the atoms are arranged in order in each branch*/
           for(jj=0;jj<=2;jj++){
              subch[i-1][j-1][k-1][jj]=subst[i].pos[atomname][jj];
           }
         }
     }/*j*/
  }/* i*/
  for(i=0;i<=2;i++) {/* add 3 Angstoms to min and max*/
      xsmax[i]=xsmax[i]+3.0;
      xsmin[i]=xsmin[i]-3.0;
  }
}    

#if 1
void myconnect(struct substrate* subst,int side[AA][NBRN][NPTS],int snum,int subname) /*snum=substrate number and subname=101*/
/*conects points on basis of distance<bondlenght; cmat[i][j] is list of pts
  j connected to atom i<j; bond() as recursion*/
{
    int i,j,natot,bi;
    float  bsqr;
    natot=subst[snum].pos[0][0]; /*total # of atoms*/
    bsqr=tsqr(bondlength);
    for(i=1;i<=natot;i++)
       cmat[i][0]=0; /*initally no connections to i*/
    side[subname][0][0]=0; /*initially no branchs*/
    for(i=1;i<=natot;i++) 
       for(j=i+1;j<=natot;j++){ /*start branch*/
         side[subname][0][1]=1; /* = initial # of points on any branch*/
         bi=0; /*indicator for new bond formed*/
         bond(subst,side,snum,subname,natot,bsqr,i,j,&bi);
         if(bi==1) /*indicates formed new branch*/
            side[subname][ side[subname][0][0] ][0]=side[subname][0][1];
       }
}

#endif
      
void bond(struct substrate* subst,int side[AA][NBRN][NPTS],int snum,int subname,int natot,float bsqr,int i,int i1,int *bi) 
//int i,i1,natot,snum,*bi;
//float bsqr;
/* checks to see if atoms i, and i1 are <bondlength, if so, it checks
   whether have already been connected, it not, it forms new bond and then
  by recursion, starts over with i1; if a bond is found at one level, then
   do not try other atoms*/
{     
      int j,k,npts,bf;
      float dist;
      dist=tsqr(subst[snum].pos[i][0]-subst[snum].pos[i1][0]) +
             tsqr(subst[snum].pos[i][1]-subst[snum].pos[i1][1])+ 
             tsqr(subst[snum].pos[i][2]-subst[snum].pos[i1][2]);
      if(dist<bsqr) { /*connect*/
           for(k=1;k<=cmat[i][0];k++) /*test whether already connected*/
              if(cmat[i][k]==i1)
                   goto precon;           
           *bi=1; /*indicates bond formed*/
           npts= side[subname][0][1]+1; /*not previous branch, make bond*/
           side[subname][0][1]=npts;
           if(npts==2){ /*first point on brach*/
               side[subname][0][0]=side[subname][0][0]+1;/*increase braanch # by 1*/
               side[subname][ side[subname][0][0] ][1]=i; /*first point*/
           }
           side[subname][ side[subname][0][0] ][npts]=i1; /*atom # on branch*/
           cmat[i][0]=cmat[i][0]+1;/*tot # pts connected to i*/ 
           cmat[i][ cmat[i][0] ]=i1; /*keeps track of atoms connected to atom i*/
           for(j=i1+1;j<=natot;j++){ /*recursion*/
              bf=0; /*indicator for new bond formed*/
              bond(subst,side,snum,subname,natot,bsqr,i1,j,&bf);
              if(bf==1) /*bond formed*/
                 goto precon; /*end branch*/
           }
     } /*dist<bsqr*/
     precon:;
}
    
void assignresnum(char nres[],char nm[],int *resn)
{
   if(strcasecmp(nm,"ALA")==0) *resn=ALA;
     else if(strcasecmp(nm,"ARG")==0) *resn=ARG;
     else if(strcasecmp(nm,"ASN")==0) *resn=ASN;
     else if(strcasecmp(nm,"ASP")==0) *resn=ASP;
     else if(strcasecmp(nm,"CYS")==0) *resn=CYS;
     else if(strcasecmp(nm,"CYH")==0) *resn=CYS; /* O uses CYH for cysteine */
     else if(strcasecmp(nm,"GLU")==0) *resn=GLU;
     else if(strcasecmp(nm,"GLN")==0) *resn=GLN;
     else if(strcasecmp(nm,"GLY")==0) *resn=GLY;
     else if(strcasecmp(nm,"HIS")==0) *resn=HIS;
     else if(strcasecmp(nm,"ILE")==0) *resn=ILE;
     else if(strcasecmp(nm,"LEU")==0) *resn=LEU;
     else if(strcasecmp(nm,"LYS")==0) *resn=LYS;
     else if(strcasecmp(nm,"MET")==0) *resn=MET;
     else if(strcasecmp(nm,"PHE")==0) *resn=PHE;
     else if(strcasecmp(nm,"PRO")==0) *resn=PRO;
     else if(strcasecmp(nm,"SER")==0) *resn=SER;
     else if(strcasecmp(nm,"THR")==0) *resn=THR;
     else if(strcasecmp(nm,"TRP")==0) *resn=TRP;
     else if(strcasecmp(nm,"TYR")==0) *resn=TYR;
     else if(strcasecmp(nm,"VAL")==0) *resn=VAL;
     else if(strcasecmp(nm,"SOL")==0) *resn=SOL;
     else if(strcasecmp(nm,"CIT")==0) *resn=CIT;
     else {
      fprintf(fpout,"%s is not in dictionary, using connect routine\n",nm);
      *resn=101; /*larger than any value in dictionary*/
     }
}

void assignatomnum(char nres[],char resname[],char na[],int *atomn)
{
   if(strcasecmp(na,"C")==0) *atomn=C;
     else if(strcasecmp(na,"CA")==0) *atomn=CA;
     else if(strcasecmp(na,"CB")==0) *atomn=CB;
     else if(strcasecmp(na,"CG")==0) *atomn=CG;
     else if(strcasecmp(na,"CD")==0) *atomn=CD;
     else if(strcasecmp(na,"CD1")==0) *atomn=CD1;
     else if(strcasecmp(na,"CD2")==0) *atomn=CD2;
     else if(strcasecmp(na,"CE")==0) *atomn=CE;
     else if(strcasecmp(na,"CE1")==0) *atomn=CE1;
     else if(strcasecmp(na,"CE2")==0) *atomn=CE2;
     else if(strcasecmp(na,"CE3")==0) *atomn=CE3;
     else if(strcasecmp(na,"CG")==0) *atomn=CG;
     else if(strcasecmp(na,"CG1")==0) *atomn=CG1;
     else if(strcasecmp(na,"CG2")==0) *atomn=CG2;
     else if(strcasecmp(na,"CH1")==0) *atomn=CH1;
     else if(strcasecmp(na,"CH2")==0) *atomn=CH2;
     else if(strcasecmp(na,"CZ")==0) *atomn=CZ;
     else if(strcasecmp(na,"CZ1")==0) *atomn=CZ1;
     else if(strcasecmp(na,"CZ2")==0) *atomn=CZ2;
     else if(strcasecmp(na,"CZ3")==0) *atomn=CZ3;
     else if(strcasecmp(na,"N")==0) *atomn=N;
     else if(strcasecmp(na,"ND1")==0) *atomn=ND1;
     else if(strcasecmp(na,"ND2")==0) *atomn=ND2;
     else if(strcasecmp(na,"NE")==0) *atomn=NE;
     else if(strcasecmp(na,"NE1")==0) *atomn=NE1;
     else if(strcasecmp(na,"NE2")==0) *atomn=NE2;
     else if(strcasecmp(na,"NH1")==0) *atomn=NH1;
     else if(strcasecmp(na,"NH2")==0) *atomn=NH2;
     else if(strcasecmp(na,"NZ")==0) *atomn=NZ;
     else if(strcasecmp(na,"O")==0) *atomn=O;
     else if(strcasecmp(na,"OD")==0) *atomn=OD;
     else if(strcasecmp(na,"OD1")==0) *atomn=OD1;
     else if(strcasecmp(na,"OD2")==0) *atomn=OD2;
     else if(strcasecmp(na,"OE1")==0) *atomn=OE1;
     else if(strcasecmp(na,"OE2")==0) *atomn=OE2;
     else if(strcasecmp(na,"OG")==0) *atomn=OG;
     else if(strcasecmp(na,"OG1")==0) *atomn=OG1;
     else if(strcasecmp(na,"OG2")==0) *atomn=OG2;
     else if(strcasecmp(na,"OH")==0) *atomn=OH;
     else if(strcasecmp(na,"SD")==0) *atomn=SD;
     else if(strcasecmp(na,"SG")==0) *atomn=SG;
     else if(strcasecmp(na,"O11")==0) *atomn=O11;
     else if(strcasecmp(na,"O12")==0) *atomn=O12;
     else if(strcasecmp(na,"C1")==0) *atomn=C1;
     else if(strcasecmp(na,"C2")==0) *atomn=C2;
     else if(strcasecmp(na,"C3")==0) *atomn=C3;
     else if(strcasecmp(na,"O3")==0) *atomn=O3;
     else if(strcasecmp(na,"C4")==0) *atomn=C4;
     else if(strcasecmp(na,"O41")==0) *atomn=O41;
     else if(strcasecmp(na,"O42")==0) *atomn=O42;
     else if(strcasecmp(na,"C5")==0) *atomn=C5;
     else if(strcasecmp(na,"C6")==0) *atomn=C6;
     else if(strcasecmp(na,"O61")==0) *atomn=O61;
     else if(strcasecmp(na,"O62")==0) *atomn=O62;
     else fprintf(stderr," ERROR IN ATOM NAME IN RESIDUE # %s, %s; no match for %s\n",nres,resname,na);
}

float tsqr(float x)
{
   return x*x;
}

void data_side(int side[AA][NBRN][NPTS],int brdraw[AA][NBRN])
/* For amino acids, the side chain data is listed as (for each branch):
            # of branchs; #of atoms in branch; ordered list of atoms in branch.
   For substrates, whole molecule is considered as series of branchs
   a[RES][list#] to side[RESNAME][..][..] where
     side[name][0][0]=# of branches in  residue name;
     side[name][i>0][0]=# of points in ith branch;
     side[name][i>0][j>0]= atom name of jth point in ith branch
  side is also used for writing pdb model, but cannot use
      duplicate atoms, thus must modify for ring a.a., use
      brdraw for writing
  
*/ 
{
  int i,j,k,count,npts,nbrn,resname; 
  for(i=0;i<SOL;i++) { /* read everything upto water, for substrates,
             whole molecule treated as sidechain */
     resname=data[i][0];
     nbrn=data[i][1];
     side[resname][0][0]=nbrn;
     count=1; 
     for(j=1;j<=nbrn;j++){
       count=count+1;      
       npts=data[i][count];
       side[resname][j][0]=npts;
       for(k=1;k<=npts;k++) {
           count=count+1;
           side[resname][j][k]=data[i][count];
       }
     }
  }
   for(i=ALA;i<=VAL;i++){ /* modifify side for writing ring aa, first copy to brdraw*/
      for(j=1;j<=side[i][0][0];j++)
          brdraw[i][j]=side[i][j][0];
   }
   brdraw[PHE][1]=8; /* modify following a.a. */
   brdraw[PRO][1]=4;
   brdraw[TRP][2]=1;
   brdraw[TRP][1]=11;
   brdraw[TYR][2]=3;
   brdraw[HIS][1]=7;
} 

void copyonefittolev(struct pepfit& fit,struct pdbfile& pdb,
      int fitres,int segn,int resn)
/*
This assumes copying to residue to allready existing (same residue type) PDB
    thus do not need resdat[][][], etc.
 copies resnum fitres in fit to lev (=pdb) segement and resn
  the restype = -fit.nside[j]  (should be >0)
NOTE: each time call - run data_side - to generate side, etc
  Faster if run just once and save 
*/
{
    int i,ii,j,jj,k,npts,nbrn,resname,count,atomname,rottype;
    int atomchnum;
    int side[AA][NBRN][NPTS];
    float pos[OHH+2][3];

    i = resn;
    resname = -fit.nside[fitres];
    nbrn=data[resname-1][1];
    side[resname][0][0]=nbrn;
    count=1; 
    for(j=1;j<=nbrn;j++){
      count=count+1;      
      npts=data[resname-1][count];
      side[resname][j][0]=npts;
      for(k=1;k<=npts;k++) {
         count=count+1;
         side[resname][j][k]=data[resname-1][count];
      }
    }// j

    for(j=0;j<=2;j++)
        pdb.mainch[segn][i][0][j]=fit.pos[fitres][N][j];
    for(j=0;j<=2;j++)
        pdb.mainch[segn][i][1][j]=fit.pos[fitres][CA][j];
    for(j=0;j<=2;j++)
       pdb.mainch[segn][i][2][j]=fit.pos[fitres][C][j];
  /*write out sidechains*/
    count = 0;
    rottype = resname;
    for(ii=1;ii<=nbrn;ii++){ // copy each side chain branch
        for(jj=0;jj<side[rottype][ii][0];jj++){
           atomchnum=side[rottype][ii][jj+1];
           if(atomchnum>5){ // not part of mainchain or Cb
              count = count+1;
              for(j=0;j<=2;j++){
                 pos[atomchnum][j]= fit.pos[fitres][count+5][j];
              }
           }
           else{ // Cb
              for(j=0;j<=2;j++)
                 pos[atomchnum][j]= fit.pos[fitres][atomchnum][j];
           }
        }// jj
    }// ii


    for(j=1;j<=nbrn;j++) { /*jth branch*/
          npts=side[resname][j][0];
          for(k=1;k<=npts;k++) {/* kth point*/
             atomname=side[resname][j][k];             
             for(jj=0;jj<=2;jj++){ /* x,y,z */
//printf("Before: pdb.sidech[%d][%d][%d][%d-1][%d] = %5.3f\n",
//   segn,i,j,k,jj,pdb.sidech[segn][i][j][k-1][jj]);
                pdb.sidech[segn][i][j][k-1][jj]=pos[atomname][jj];
//printf("After: pdb.sidech[%d][%d][%d][%d-1][%d] = %5.3f\n",
//   segn,i,j,k,jj,pdb.sidech[segn][i][j][k-1][jj]);
             }/* jj */
          }/* k */
    }/* j */
    /* add C-O branch*/
    for(jj=0;jj<=2;jj++) 
        pdb.sidech[segn][i][j][0][jj]=fit.pos[fitres][C][jj];
    for(jj=0;jj<=2;jj++) 
        pdb.sidech[segn][i][j][1][jj]=fit.pos[fitres][O][jj];
}


void copylevtofit(struct pepfit& fit,struct pdbfile& pdb,char pdbfile[])
/* copies entire lev to fit
*/
{
   int segn;

   for(segn=0;segn<pdb.resdat[0][0][0];segn++){ // all segments
       copyreglevtofit(fit,pdb,segn,1,pdb.resdat[segn+1][0][0]);
   }
}  

void writelevtopdb(struct pdbfile& pdb, const char pdbnew[])
/* Routine to write the screen pdb (in lev format) to
    file pdbnew.
  REads through original PDB file (= pdb.name) use to make lev,
   for each residu: (segn,resn)
   1)  Checks if it has been changed - USED AGAIN (15JUN00) NOT USED ANYMORE (18may00)
          indicated by resdat[segn+1][resn][4]
   2) If changed, reads and writes main chain exactly.
        for Side chain atoms, reads each line, change only 
           the pos, and writes again.
*/
{
   int j,jj,k,atomnum,count;
   int nres,nseg,rnum,nbrn,npts;
   int inew,itotal;

   FILE *fp,*fpout;
   char lresnumch[NCHAR],oldresnumch[NCHAR]; 
   char dataline[200],atomch[30],atomnumch[NCHAR],atomname[NCHAR];
   char resname[NCHAR],seg[NCHAR],oldseg[NCHAR];
   float sigma,bfact,x[3];
   float pos[OHH+2][3];
   int readatom;

   strcpy(oldresnumch,"");
   nseg= -1;
   nres=0;
   strcpy(oldseg,"");
   strcpy(seg,"");/* need this because if seg part of data line missing, will use whatever is intial value*/
  if( (fp=fopen(pdb.name,"r") ) ==NULL) {
      fprintf(stderr,"error opening file %s\n",pdb.name);
      exit(0);
  }
  fpout = fopen(pdbnew,"w");
  inew = itotal = 0;
   readatom = 0;
  while(  fgets(dataline,199,fp) != NULL) {
    if(strncmp("TER",dataline,3)==0){ // If line starts with TER, new segment
       readatom = 0;
       fprintf(fpout,"%s",dataline); // Print TER line
    }
    else if(strncmp("ATOM",dataline,4)!=0){ // If does not begin with ATOM, just write again
       fprintf(fpout,"%s",dataline); 
    }
    else { // begins with ATOM
        if(readatom==0){
           readatom = 1;
           nseg = nseg+1;
           nres = 0;
        }
        if(dataline[21]!=' ') /* then standard pdb format */
            sscanf(dataline,"%s%s%s%s%s%s%8f%8f%8f",atomch,atomnumch,atomname,
               resname,seg,lresnumch, &x[0],&x[1],&x[2]);
        else{ /* xplor format */
           sscanf(dataline,"%s%s%s%s%s%8f%8f%8f%6f%6f%s",atomch,atomnumch,atomname,resname,lresnumch,
                                     &x[0],&x[1],&x[2],&sigma,&bfact,seg);
        }
#if 0
       if( strcmp(seg,"") == 0) /* segname blank*/
           nseg=0;
       else { /*sement name not blank*/
            if( (seg[0]<'A')||(seg[0]>'Z') ) /* 65 is character A, if not A...Z then assume null */
                 nseg=0;
            else if( strcmp(seg,oldseg)!=0) {/*change in segment name,*/
               strcpy(oldseg,seg);
               nseg=nseg+1;
               nres=0;        
            }
        } 
#endif
       /*discard H atoms and drop OT2,change OT1 to O (xplor format
             in pdb format, drop OT, and discard waters*/
       if( (strncmp(atomname,"H",1)!=0)&& (strcmp(atomname,"OT")!=0)&&
           (strcmp(atomname,"OT2") !=0)&&(strcmp(atomname,"OXT") !=0)&&
           (strncmp(atomname,"1H",2)!=0)&&(strncmp(atomname,"2H",2)!=0)&&
            (strncmp(atomname,"3H",2)!=0)  &&
            (strcmp(resname,"SOL")!=0)&&(strcmp(resname,"WAT")!=0) ){
// then read line
           if(strcmp(atomname,"OT1") ==0)
                  strcpy(atomname,"O");
           if(strcmp(resname,"ILE")==0) /* in ILE, allow either CD or CD1*/
              if(strcmp(atomname,"CD")==0)
                 strcpy(atomname,"CD1");
           if(strcmp(oldresnumch,lresnumch)!=0){ /*new residue - check distance of last 2 residues*/
//printf("ires = %d  ifit = %d  fit[ifit].pos[ires][1][0] = %5.3f\n",ires,ifit, fit[ifit].pos[ires][1][0]);
              strcpy(oldresnumch,lresnumch);
              assignresnum(lresnumch,resname,&rnum);/*if not in dictionary,assigned rnum=101*/
              if(rnum<=VAL){ // then add new residue
                nres=nres+1;
              }
#if 1 //15jun00 - use 4 in place of 100 and increase NBRN from 4 to 5 Not used anymore (18may00)  
              itotal = itotal+1;
              if(pdb.resdat[nseg+1][nres][4] == 1){ // Then has been changed  -Read all side atoms 
                inew = inew+1;             
                nbrn=data[rnum-1][1]; 
                side[rnum][0][0]=nbrn;
                count=1;
                for(j=1;j<=nbrn;j++) { /*jth branch*/
                   count=count+1;      
                   npts=data[rnum-1][count];
                   for(k=1;k<=npts;k++) {/* kth point*/
                     count=count+1;
                     atomnum=data[rnum-1][count];
                     for(jj=0;jj<=2;jj++){ /* x,y,z */
                        pos[atomnum][jj]=pdb.sidech[nseg][nres][j][k-1][jj];
                     }//jj
                   } // k
                 }//j
              }// has been changed
#endif
            } // new res
#if 1 //15jun00 - use 4 in place of 100 and increase NBRN from 4 to 5 Not used anymore (18may00) Always = 0
           if(pdb.resdat[nseg+1][nres][4] == 1){ // Then has been changed  -Read all side atoms                         
              assignatomnum(lresnumch,resname,atomname,&atomnum);               
              checkatomname(side,lresnumch,rnum,&atomnum); //checks names and  make correct to ILE
              if(atomnum<=4){ // not a side chain, just write
                   fprintf(fpout,"%s",dataline);
              }
              else{ // a side chain atom - write out with new  positions
                 fprintf(fpout,"ATOM  %5s  %-3s %3s  %4s     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
                     atomnumch,atomname,resname,lresnumch,
                     pos[atomnum][0],pos[atomnum][1],pos[atomnum][2],sigma,bfact,seg);
              }
           } // has been changed
           else{ // has not been changed, just print dataline
                   fprintf(fpout,"%s",dataline);
           }
#endif
//           fprintf(fpout,"%s",dataline); //Use this line if comment out above and write evry line
       } /* if not H or water */ 
       else{ // just print line
          fprintf(fpout,"%s",dataline);
       } 
      } /* if ATOM*/
   } // while
   printf(" FRACTION OF RESIDUES IMPROVED = %5.3f\n",(float)inew/itotal);
   fclose(fp);
   fclose(fpout);
}
 

void copyreglevtofit(struct pepfit& fit,struct pdbfile& pdb,
      int segn,int lev0,int lev1)
/*
  Makes a fit from region lev0 to lev1  of a lev (PDB) file.
     segn is the segment number of PDB.
NOTE: data[resname][] starts at zero instead of 1 
  thus must shift resname down by one
For fit from 1 to N need C and O  from fit[0]
Thus copies to fit 0,1,2,.. fitres = fit,nres-1,...
 the set fit.nres = fit.nres-1 at end
*/
{
    int i,ii,j,jj,k,npts,nbrn,resname,count,fitres,atomname;
    int side[AA][NBRN][NPTS];
    float pos[OHH+2][3];
    char dummy[NCHAR];

    ii = segn;
    fit.nres = 0;
    for(i = lev0;i<=lev1;i++){
     if( (i<1)||(i>pdb.resdat[segn+1][0][0]) ){
       fit.nres = fit.nres+1;
       fitres = fit.nres -1;
       resname = GLY;
       for(jj=1;jj<=4;jj++) // mainchain atoms
         for(k=0;k<=2;k++)
          fit.pos[fitres][jj][k] = 0.0;
       fit.itype[fitres] = 0;
       fit.nside[fitres]= resname;
       fit.nres = fitres+1;
     }
     else{
       fit.nres = fit.nres+1;
       fitres = fit.nres -1;
       assignresnum(dummy,pdb.rname[ii][i],&resname);
//       fit.nside[fitres] = -resname;  // done in copypostofit 
       nbrn=data[resname-1][1];
       side[resname][0][0]=nbrn;
       count=1; 
       for(j=1;j<=nbrn;j++){
          count=count+1;      
          npts=data[resname-1][count];
          side[resname][j][0]=npts;
          for(k=1;k<=npts;k++) {
            count=count+1;
            side[resname][j][k]=data[resname-1][count];
          }
       }// j
 
       for(j=0;j<=2;j++)
          pos[N][j]=pdb.mainch[ii][i][0][j];
       for(j=0;j<=2;j++)
          pos[CA][j]=pdb.mainch[ii][i][1][j];
       for(j=0;j<=2;j++)
          pos[C][j]=pdb.mainch[ii][i][2][j];
       /*write out sidechains*/
       for(j=1;j<=nbrn;j++) { /*jth branch*/
          npts=side[resname][j][0];
          for(k=1;k<=npts;k++) {/* kth point*/
             atomname=side[resname][j][k];             
             for(jj=0;jj<=2;jj++){ /* x,y,z */
                pos[atomname][jj]=pdb.sidech[ii][i][j][k-1][jj];
            }/* jj */
          }/* k */
        }/* j */
        /* add C-O branch*/
//        for(jj=0;jj<=2;jj++) 
//           pos[C][jj]=pdb.sidech[ii][i][j][0][jj];
        for(jj=0;jj<=2;jj++) 
           pos[O][jj]=pdb.sidech[ii][i][j][1][jj];
        copypostofit(fit,pos,fitres,resname,dummy);
    }
   }/*for i*/
   fit.nres = fit.nres -1;
}
