/* fittopdb.c++  
  Routines for converting fit file to pdb and pdb to fit
Modified from thinfittopdb.c in old maid with following major changes:
I) glmaid fit has major change from maid fit
  glmaid fit: fit.pos[residue #][atom #][x,y,z
     for(j=0;j<fit[i].nres;j++){ // each residue
       for(ii=1;ii<=4+fit[i].nside[j];ii++) // all atom positions per reisude 
          for(k=0;k<3;k++) // x y z 
            fprintf(fp," %f",fit[i].pos[j][ii][k]);  
Fit 1   startendpt 281 4177 
 9
 1 1 1 1 1 1 1 1 1
 0 5 1 3 4 2 5 2 0
 19.662111 28.747742 -5.918052 19.074137 30.149807 -5.775826 18.406065 30.656706 -6.710850 0.000000 0.000000 0.000000
 18.845217 32.098248 -4.304376 19.742054 33.132843 -4.980409 19.346928 34.315643 -5.130068 19.325460 30.748093 -4.627499 18.882374 32.303238 -2.791654 18.480314 33.722416 -2.396669 17.033924 33.841438 -2.169179 16.647137 35.206711 -1.789197 17.795862 35.994259 -1.322117
 21.917305 33.526268 -6.042242 21.463106 33.818119 -7.470617 21.744770 34.914051 -8.015720 20.917912 32.682606 -5.373356 23.259171 32.798157 -6.073557
 20.256153 32.944775 -9.418571 19.208254 34.053215 -9.489278 19.162239 34.822002 -10.481438 20.777681 32.850643 -8.048494 19.617672 31.614824 -9.812655 19.744324 31.340788 -11.309516 18.457964 31.479725 -12.005169
 17.322203 35.118343 -8.339370 17.930759 36.518650 -8.316484 17.368126 37.466095 -8.919248 18.394117 34.120403 -8.453701 16.541725 34.881382 -7.048466 15.569391 33.710743 -7.174578 15.138175 33.209766 -5.862696 14.966039 31.750887 -5.860184
 19.787104 37.894947 -7.495646 20.291208 38.342342 -8.865844 20.069138 39.508896 -9.274956 19.053375 36.629745 -7.633000 20.974630 37.695316 -6.556683 21.424585 39.006680 -5.916684
 21.513863 37.666580 -10.881392 20.402750 38.116825 -11.827141 20.672485 38.804996 -12.842624 20.952118 37.425377 -9.545645 22.141310 36.378654 -11.409865 22.865990 36.597641 -12.735992 23.469604 35.358635 -13.244393 23.224686 34.219685 -12.349538 24.281866 34.075962 -11.339749
 17.991606 38.064091 -12.271673 17.437878 39.408916 -11.806317 17.141760 40.298038 -12.642567 19.188297 37.729668 -11.488054 16.933487 36.981384 -12.072039 15.534399 37.475407 -12.432903
 16.792234 40.757523 -9.862511 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 17.307358 39.538528 -10.500207

NOTE from above fit file:
  First residue - corresponds to Ca (nterm)-C-O (no side chain)
  Middle reisdue  N-Ca-C-O
  Last residue - N-Ca (cterm) no side chain

  maid fit: pos[atom #][residue #][x,y,z]
   for(i=nend[0];i<=nend[1];i++){
     ipos=ipos+1;
     for(j=1;j<=6;j++)
       for(k=0;k<=2;k++){
          fprintf(fp," %f %f",pos[j][ipos][k],vel[j][ipos][k]);

II) In glmaid - number of side chain atoms varies from 0 to 5,
     when make pdb - if >2 trucate at 2 and write as serine
       if 1 side chain - write as alanine
       if 0 side chain - write as glycine

III) Handle segments differently - place a "TER Residue #" after each 
       fit fragment.

IV) Only option - write entire fit  (do not have option of choosing one fit to write)
V)  Only option - write only complete residues:
     This means drop Ca, C and O and nterm end and N and Ca at Cterm end
VI) getyn is not used in glmaid
VII) change definitino of aaname:
 in maid:
   strcpy(aaname[0],"END");
   strcpy(aaname[1],"SER");
   strcpy(aaname[2],"GLY");
 in glmaid:
   strcpy(aaname[0],"GLY");
   strcpy(aaname[1],"ALA");
   strcpy(aaname[2],"SER");
VIII) in maid - each fit is a separate file, in maid - one file
*/
#include "peptide.h"
#include "drawwindow.h"  // also calls global.h
#include "trace.h"
#include "pdbtolev.h"
#include "math.h"
#include "fittopdb.h"
#include "string.h"
#include "stdlib.h" 

extern char*  nameres[30];
extern char* nameatom[100];
extern int side[AA][NBRN][NPTS];  // defined globally (originally in pdbtolev.c++)  because needed in assignseg.c++ 


void copypostofit(struct pepfit& fit,float pos[OHH+2][3],int ires,int rnum,char seg[])
{
    int i,ii,j,k,count,atomnum;

#if 0 // Now - all fits pre- allocated
    if(ires >=fit.nalloc)
       nalloc = ires+20;
    else
       nalloc = fit.nalloc;
    allocatefit(fit,nalloc);
#endif
    fit.nside[ires] = -rnum;
    fit.itype[ires] = 0;// default - reassigned below
    fit.nres = ires+1;
    for(k=0;k<=2;k++){
      for(i=1;i<=4;i++) // mainchain atoms
        fit.pos[ires][i][k] = pos[i][k];
      count = 0;
      for(ii=1;ii<=side[rnum][0][0]+1;ii++){ // copy each side chain branch
          for(j=0;j<side[rnum][ii][0];j++){
             atomnum=side[rnum][ii][j+1];
             if(atomnum>5){ // not part of mainchain or Cb
               count = count+1;
               fit.pos[ires][count+5][k]=pos[atomnum][k];
             }
             else{// = 5 = Cb
                fit.pos[ires][atomnum][k]=pos[atomnum][k];
             }
          } // j
      } // ii
   }//k
   if(strcmp(seg,"NASS")==0){ // if not assigned, set val1 and side type
     fit.val1 = 0;
     if(rnum==GLY)
       fit.nside[ires] = 0;
     else if(rnum == ALA)
       fit.nside[ires] = 1;
     else if(rnum == SER)
       fit.nside[ires] = 2;
   }
}


void sizefit(struct geometry& geo,char fitpdbfile[],int size[NFIT])
/* 
  Similar to writepdbtofit, but just goes through entire pdb file,
    finds gaps and returns the size of each fit.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int n,k,ifit,ires,rnum,atomnum;
   char lresnumch[NCHAR],oldresnumch[NCHAR]; 
   char dataline[200],atomch[30],atomnumch[NCHAR],atomname[NCHAR],
       resname[NCHAR],seg[NCHAR];
   float x[3],sigma,bfact,dist,pos[OHH+2][3],oldpos[OHH+2][3];
   FILE *fp;
  
   float gapdist = 3.0;  // max allowed distance before assume a "gap"
   if( (fp=fopen(fitpdbfile,"r") ) ==NULL) {
      fprintf(stderr,"error opening file %s\n",fitpdbfile);
      exit(0);
   }
   ifit=1; // first fit 
   ires = 0; // first residue number 
   strcpy(oldresnumch,""); // set blank so will no when new name read
   n=0;
   while(  fgets(dataline,199,fp) != NULL) {
    if(strncmp("ATOM",dataline,4)==0){ /*Read only if begins with ATOM*/
        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);
       }
       /*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*/
             strcpy(oldresnumch,lresnumch);
             n=n+1;
             if(n==2){
                for(k=0;k<=2;k++){
                  oldpos[2][k] = pos[2][k];
                  oldpos[4][k] = pos[4][k];
                }
             }
             if(n>2){ // check distance
                dist = distv(oldpos[2],pos[4]); // C - N dist
                if(dist>gapdist){// then pos is start of another fit
printf("A gap!!!!\n");
                   size[ifit] = ires +1;
                   ifit = ifit+1;
                   ires =  0;
                }
                else  // add pos to prev. fit and classify angle type
                   ires =  ires+1;
                for(k=0;k<=2;k++){
                      oldpos[2][k] = pos[2][k];
                      oldpos[4][k] = pos[4][k];
                }
            }
            assignresnum(lresnumch,resname,&rnum);/*if not in dictionary,assigned rnum=101*/
           }// new residue
           if(rnum<=VAL){ // then add new residue
               assignatomnum(lresnumch,resname,atomname,&atomnum);
               if( (atomnum ==2)||(atomnum==4) ){ // save N and C positin              
                  checkatomname(side,lresnumch,rnum,&atomnum); //checks names and  make correct to ILE
                  for(k=0;k<=2;k++)
                      pos[atomnum][k]= x[k];
               }
           }
       } /* if not H or water */ 
      } /* if ATOM*/
   } // while
   if(oldpos[2][0]!=pos[2][0]){ // then read last residue
        dist = distv(oldpos[2],pos[4]); // C - N dist
        if(dist>gapdist){// then pos is start of another fit
printf("A gap!!!!\n");
            size[ifit] = ires +1;
            ifit = ifit+1;
            size[ifit] = 2; // last residue is by itself 
        }
        else{  // add pos to prev. fit
           ires= ires+1;
           size[ifit] = ires+1;
           size[0] = ifit; // total number of fits
       }
   }
   fclose(fp);
   for(int i=1;i<=size[0];i++)
     printf("For fit %d size = %d\n",i,size[i]);
}


void writepdbtofit(struct geometry& geo,struct pepfit fit[NFIT],int fitnum[NFIT],
            char fitpdbfile[])
/* 
NEW: assignes val1 = starting sequence number 
Creates set of fits from pdbfile and write to standard fitfile
  Reads pdbfile and creates fit, one res at a time.
   Checks the C[i] to N[i+1] distance and if too long, assumes break
     ends one fit and starts another
  In order to find C and N must read and store last 2 residues
  Read pdb into fit as if no gap, after read, check distance, if gap
    then copy this last residue into first residue of new fit.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,n,k,ifit,ires,rnum,atomnum;
   int origntype,size[NFIT]; // Goes through entire file, finds size of each fit
   char lresnumch[NCHAR],oldresnumch[NCHAR]; 
   char dataline[200],atomch[30],atomnumch[NCHAR],atomname[NCHAR],
       resname[NCHAR],seg[NCHAR],segchar[NCHAR];
   float x[3],sigma,bfact,dist,pos[OHH+2][3],phi,psi;
   FILE *fp;
  
   float gapdist = 3.0;  // max allowed distance before assume a "gap"
   origntype = geo.ntypes;
   geo.ntypes = 7; // Include gammar and gammal
   sizefit(geo,fitpdbfile,size);
   for(i=1;i<=size[0];i++){
     if(i> NFIT-1){
       printf("ERROR:  i > NFIT in writepdbtofit\n");
       goto quitrout;
     }
     allocatefit(fit[i],size[i]);
   }     
   if( (fp=fopen(fitpdbfile,"r") ) ==NULL) {
      fprintf(stderr,"error opening file %s\n",fitpdbfile);
      exit(0);
   }
   ifit=fitnum[0]=1; // first fit 
   ires = 0; // first residue number 
   strcpy(oldresnumch,""); // set blank so will no when new name read
   n=0;
   while(  fgets(dataline,199,fp) != NULL) {
    if(strncmp("ATOM",dataline,4)==0){ /*Read only if begins with ATOM*/
        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);
       }
       /*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);
             n=n+1;
             if(n==2){
                fit[ifit].val1 = -atoi(lresnumch);
                fit[ifit].val2 = 0;
                copypostofit(fit[ifit],pos,ires,rnum,seg);
                fitnum[ifit] = ifit;
             }
             if(n>2){ // check distance
                dist = distv(fit[ifit].pos[ires][2],pos[4]); // C - N dist
                if(dist>gapdist){// then pos is start of another fit
printf("A gap!!!!\n");
                   ifit = fitnum[0]= fitnum[0]+1;
                   fit[ifit].val1 = -atoi(lresnumch);                   
                   fit[ifit].val2 = 0;
                   ires =  0;
                }
                else  // add pos to prev. fit and classify angle type
                   ires =  ires+1;
                copypostofit(fit[ifit],pos,ires,rnum,segchar); // segchar is from prevoius line
                fitnum[ifit] = ifit;
            }
             assignresnum(lresnumch,resname,&rnum);/*if not in dictionary,assigned rnum=101*/
           }// new residue
           if(rnum<=VAL){ // then add new residue
               assignatomnum(lresnumch,resname,atomname,&atomnum);               
               checkatomname(side,lresnumch,rnum,&atomnum); //checks names and  make correct to ILE
               for(k=0;k<=2;k++)
                   pos[atomnum][k]= x[k];
               strcpy(segchar,seg);
           }
       } /* if not H or water */ 
      } /* if ATOM*/
   } // while
   if(fit[ifit].pos[ires][1][0]!=pos[1][0]){ // then read last residue
        dist = distv(fit[ifit].pos[ires][2],pos[4]); // C - N dist
        if(dist>gapdist){// then pos is start of another fit
printf("A gap!!!!\n");
            ifit = fitnum[0]= fitnum[0]+1;
            ires = fit[ifit].nres = 0;
        }
        else  // add pos to prev. fit
           ires = fit[ifit].nres = ires+1;
        copypostofit(fit[ifit],pos,ires,rnum,seg);
        fitnum[ifit] = ifit;
   }
   fclose(fp);
   for(ifit = 1;ifit<=fitnum[0];ifit++){
     for(ires=1;ires<fit[ifit].nres-1;ires++){ // then classify phi/psi angles for ires
          findphipsifit(fit[ifit],ires-1,&phi,&psi);
          nclassifyphipsi(geo,phi,psi,&fit[ifit].itype[ires]); // Uses geo.sheet parameters
//printf("phi = %5.3f psi = %5.3f type = %d\n",phi*180/PI,psi*180/PI,fit[ifit].itype[ires]);
     }
   }              
   savefit();
   geo.ntypes = origntype;
   quitrout:;
}


void writeonefittopdb(FILE *fp,struct griddata& den1,struct geometry& geo,
            struct pepfit fit[],int fitnum[],
            char fitpdbfile[],int ifit,char reschar[],int *resnumin,int *atomnumin,
            char segchar[],int writeter)
/*
If writeter = 1, then write TER at end of fit
New - 27jun00 - add glycine at each end
   Use same type (sheet or helix) as fit, and just use ideal phi/psi angle.
Write fit # ifit to fitpdbfile, starting with resnum+1,atomnum+1
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int k,ii,j,jj,rottype,count,atomchnum,writecb,writecz,writecg;
   char atomname[10][10],aaname[10][10];
   char atomchar[10];
   float occup,tfact;
   int resnum,atomnum;

   strcpy(atomname[1],"CA");/*assign names to atoms 1 - 6 */
   strcpy(atomname[2],"C");
   strcpy(atomname[3],"O");
   strcpy(atomname[4],"N");
   strcpy(atomname[5],"CB");
   strcpy(atomname[6],"OG");
   strcpy(aaname[0],"GLY");
   strcpy(aaname[1],"ALA");
   strcpy(aaname[2],"SER");

//   strcpy(segchar,""); // default segment is blank
   occup = 1.0; // default occupancy
   tfact = 20.0; // default tempfact
   resnum = *resnumin-1; // MUST SUBTRACT 1 BECAUSE ADD GLYCINE AT NTERM END
   atomnum = *atomnumin;

     for(j=0;j<=fit[ifit].nres-1;j++){ // each complete residue - NOTE - does not include first/last residue
       if(j==0){ // add glycine geo residue at nterm end
          copyregfittogeo(geo,fit[ifit],0,4); 
          addnterm(geo,den1,fit[ifit].itype[1]); // Use same contraint as res #1
          fit[ifit].nside[j] = -GLY;
          for(k=0;k<=2;k++)// Use geo N pos made by addnterm
              fit[ifit].pos[j][4][k] = geo.pos[geo.nend[0]+1][4][k];
       } // j=0
       else if(j==fit[ifit].nres-1){ // add geo residue at cterm end
          copyregfittogeo(geo,fit[ifit],j-4,j); //NEW
          addcterm(geo,den1,fit[ifit].itype[j-1]); // Use same contraint as res #1
          fit[ifit].nside[j] = -GLY;
          for(k=0;k<=2;k++){// Use geo C and O pos made by addcterm
              fit[ifit].pos[j][2][k] = geo.pos[geo.nend[1]-1][2][k];
              fit[ifit].pos[j][3][k] = geo.pos[geo.nend[1]-1][3][k];
          }
       } 
       resnum = resnum+1;  // Start writing mainchain atoms N,Ca, C, O
       if(fit[ifit].nside[j]>=0){ // then side is just series or C atoms
           if(fit[ifit].nside[j]>2) // maximum of 2 side chain atoms - ON NEXT FIT RESIDUE
               fit[ifit].nside[j]=2;
           strcpy(reschar,aaname[ fit[ifit].nside[j] ]);
       }
       else{
           rottype = -fit[ifit].nside[j];
           strcpy(reschar,nameres[rottype]);
       }
       if(fit[ifit].nside[j]>2) // maximum of 2 side chain atoms - ON NEXT FIT RESIDUE
            fit[ifit].nside[j]=2;
       ii=4;  // N
       atomnum =atomnum+1;
       fprintf(fp,"ATOM  %5d  %-3s %3s  %4d     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
              atomnum,atomname[ii],reschar,resnum,fit[ifit].pos[j][ii][0],
              fit[ifit].pos[j][ii][1],fit[ifit].pos[j][ii][2],occup,tfact,segchar);
       ii=1;  // Ca
       atomnum =atomnum+1;
       fprintf(fp,"ATOM  %5d  %-3s %3s  %4d     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
              atomnum,atomname[ii],reschar,resnum,fit[ifit].pos[j][ii][0],
              fit[ifit].pos[j][ii][1],fit[ifit].pos[j][ii][2],occup,tfact,segchar);
       ii=2;  // C 
       atomnum =atomnum+1;
       fprintf(fp,"ATOM  %5d  %-3s %3s  %4d     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
              atomnum,atomname[ii],reschar,resnum,fit[ifit].pos[j][ii][0],
              fit[ifit].pos[j][ii][1],fit[ifit].pos[j][ii][2],occup,tfact,segchar);
       ii=3;  // O 
       atomnum =atomnum+1;
       fprintf(fp,"ATOM  %5d  %-3s %3s  %4d     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
              atomnum,atomname[ii],reschar,resnum,fit[ifit].pos[j][ii][0],
              fit[ifit].pos[j][ii][1],fit[ifit].pos[j][ii][2],occup,tfact,segchar);         
       if(fit[ifit].nside[j]>=0){ // then side is just series or C atoms
         for(jj= 5;jj<= 4+fit[ifit].nside[j];jj++){ // side chains from next residue
          atomnum =atomnum+1;
          fprintf(fp,"ATOM  %5d  %-3s %3s  %4d     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
              atomnum,atomname[jj],reschar,resnum,fit[ifit].pos[j][jj][0],
                  fit[ifit].pos[j][jj][1],fit[ifit].pos[j][jj][2],occup,tfact,segchar);
         }
       }// iside>=0
       else{ // side = assigned residue
         count = 0;
         writecb=writecg=writecz=0; //make sure only write cb once
         for(ii=1;ii<=side[rottype][0][0]+1;ii++){ // copy each side chain branch
            if(rottype==TRP && ii==2)// For arg, last branch just connects 2 previously drawn atoms
               goto nextjj;
            for(jj=0;jj<side[rottype][ii][0];jj++){
              atomchnum=side[rottype][ii][jj+1];
              strcpy(atomchar,nameatom[atomchnum]);
              if(atomchnum>5){ // not part of mainchain or Cb
                count = count+1;
                if (ii>1 && jj==0) //New 13aug00 - do not write branch atom twice
                   goto nextjj;
                if (atomchnum==CG){ // do not write cg twice
                  if( writecg==1)
                      goto nextjj;
                  else
                     writecg=1;
                }                   
                if (atomchnum==CZ){
                  if( writecz==1)
                      goto nextjj;
                  else
                     writecz=1;
                }                   
                fprintf(fp,"ATOM  %5d  %-3s %3s  %4d     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
                  atomnum,atomchar,reschar,resnum,fit[ifit].pos[j][count+5][0],
                  fit[ifit].pos[j][count+5][1],fit[ifit].pos[j][count+5][2],occup,tfact,segchar);
              }
              else if(atomchnum==CB){ // Cb
               if (writecb==0){// only write cb if have not written before
                  writecb=1;
                  fprintf(fp,"ATOM  %5d  %-3s %3s  %4d     %7.3f %7.3f %7.3f %5.2f %5.2f      %-4s\n",
                    atomnum,atomchar,reschar,resnum,fit[ifit].pos[j][atomchnum][0],
                    fit[ifit].pos[j][atomchnum][1],fit[ifit].pos[j][atomchnum][2],occup,tfact,segchar);
               }//writecb=0
              } 
              nextjj:;
            }// jj
          }// ii
        } // side = assigned residue
     } // j = res #
     if(writeter==1){// write TER at end of fit
        atomnum =atomnum+1;  // give a new atom number to TER    
        fprintf(fp,"TER  %6d  %-3s %3s  %4d\n",
              atomnum,"",reschar,resnum);
     }
     *resnumin = resnum;
     *atomnumin = atomnum;
}  

#if 1
void writefittopdb(struct griddata& den1,struct geometry& geo,struct pepfit fit[],int fitnum[],
            char fitpdbfile[])
/*
Called only if first fit doe not have segment assigned - writes all fits as
   separate segment
Assumes fit is allready read
If side chain residues have been assigned (indicated by
     fit.nside[j]<0 = - rottype)
  Then write standard pdb
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
    4+fit[i].nside[j] is side chain atom number
*/
{
   FILE *fp;
   int i,ifit,resnum,atomnum,writeter;
   char segchar[5],reschar[20];

printf("Converting  %d fits to PDB format and write to file %s\n",fitnum[0],fitpdbfile);
   atomnum = resnum = 0;
   writeter = 1; // write TER after each fit
   fp= fopen(fitpdbfile,"w");
   for(i=1;i<=fitnum[0];i++){ // write each fit fragment
     ifit = fitnum[i];
     strcpy(segchar,"NASS"); // indicates sequence not yet assigned
     writeonefittopdb(fp,den1,geo,fit,fitnum,fitpdbfile,ifit,reschar,
            &resnum,&atomnum,segchar,writeter);// return resnum = last residue written
   } // i = fit #
   fprintf(fp,"END\n");
   fclose(fp);  
}
#endif

  
#if 0  // NEW version
void writefittopdbassign(struct griddata& den1,struct geometry& geo,
            struct pepfit fit[NFIT],int fitnum[NFIT],
            char fitpdbfile[])
/*
NOTE: writeonefittopdb TO PDB ADDS A GLYCINE AT EACH END.
Assumes fit is allready read.
Sequence assignment of fit.pos[1] is =  -fit.val1.  If fit.val1 is <0
     the sequence as been assigned.
If fit.val2< -900, then segment assignment as been maid
     segment 0 = -1000;  segment 1 = -1001, etc.
If side chain residues have been assigned (indicated by
     fit.nside[j]<0 = - rottype)   Then write standard pdb
First goes through fits and looks for fits that have both sequence
  and segment assigned:
     For  each segment and write output pdb in order of a.a. number 
  A TER statement is put at end of each segment.
Then writes fits with segments not assigned.
  A TER statemente is put at end of each fit.
Then write fits with sequence not assigned.
  A TER statemente is put at end of each fit.
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
    4+fit[i].nside[j] is side chain atom number
*/
{
   FILE *fp;
   int i,ifit,resnum,atomnum,iseg,writeter,nseg[200],orderseg[200];
   char segchar[5],reschar[20];

printf("Converting  %d fits to PDB format and write to file %s\n",fitnum[0],fitpdbfile);
   writeter = 1;  // indicates write TER at end of fit
   fp= fopen(fitpdbfile,"w");
   nseg[0]=0; // counts segment numbers previous used
   for(ifit=1;ifit<=fitnum[0];ifit++){ // find all segment numbers
     if(fit[ifit].val2 < -900) {
         iseg =  -fit[ifit].val2 - 900;
         for(i=1;i<=nseg[0];i++){
            if(nseg[i]==iseg) // already used
              goto nextifit;
         }
         nseg[0]=nseg[0]+1;
         nseg[ nseg[0] ] = iseg;
     }// < -900
     nextifit:;
   }//ifit 
   orderlist(nseg,orderseg); //do not need this since already ordered?? segment numbers in order of lowest to highest   
   for(i=1;i<=orderseg[0];i++){// write segments in order
     for(ifit=1;ifit<=fitnum[0];ifit++){ // write fits sequence  and segment assignment
       if( fit[ifit].val2 < -900 ){ // segment assigned
          iseg =  -fit[ifit].val2 - 900;
          if(iseg == orderseg[i]){
             itoa(iseg,segchar);
             atomnum = 0;
             if(fit[ifit].val1<0)
                resnum = -fit[ifit].val1 -1;//Need to subtact 1 because writeonefit start at resnum+1
             else
                resnum = 0;
             writeonefittopdb(fp,den1,geo,fit,fitnum,fitpdbfile,ifit,reschar,
                 &resnum,&atomnum,segchar,writeter);// return resnum = last residue written
          }
       }
     }//ifit
   }//i
   for(ifit=1;ifit<=fitnum[0];ifit++){ // write fits sequence assignment and no segment assigned
       if( (fit[ifit].val1<0) && (fit[ifit].val2 > -900) ){ // sequence and segment assigned
          strcpy(segchar,"NSEG"); // indicates segment not yet assigned
          writeter = 1;  // indicates write TER at end of fit
          atomnum = 0;
          resnum = -fit[ifit].val1 -1;//Need to subtact 1 because writeonefit start at resnum+1
          writeonefittopdb(fp,den1,geo,fit,fitnum,fitpdbfile,ifit,reschar,
               &resnum,&atomnum,segchar,writeter);// return resnum = last residue written
       }
   }
   for(ifit=1;ifit<=fitnum[0];ifit++){ // write fits with no sequence assignment
      if( (fit[ifit].val1 >= 0)&&(fit[ifit].val2 > -900) ){
         strcpy(segchar,"NASS"); // indicates sequence and segment not yet assigned
         writeter = 1;  // indicates write TER at end of fit
         atomnum = 0;
         resnum = 0;
         writeonefittopdb(fp,den1,geo,fit,fitnum,fitpdbfile,ifit,reschar,
            &resnum,&atomnum,segchar,writeter);// return resnum = last residue written
      }
   }
   quitrout:;
   fclose(fp);  
}
#endif

#if 1  // NEWest version - 27nov00 -
void writefittopdbassign(struct griddata& den1,struct geometry& geo,
            struct pepfit fit[NFIT],int fitnum[NFIT],
            char fitpdbfile[])
/*
CALLED ONLY IF FIRST FIT HAS fit[ifit].val2 < -900 (that is segment assigned) -
   assumes all fits have segment assigned.
NEw -- Assumes fitfile already correctly ordered - with:
   segment ordered, and within segment, redidues numbered correctly
   All fits have segment assigned.
   Set writeter = 1 only for the last fit in a segment (indicates end of segment4)
NOTE: writeonefittopdb TO PDB ADDS A GLYCINE AT EACH END.
Assumes fit is allready read.
Sequence assignment of fit.pos[1] is =  -fit.val1.  If fit.val1 is <0
     the sequence as been assigned.
If fit.val2< -900, then segment assignment as been maid
     segment 0 = -1000;  segment 1 = -1001, etc.
If side chain residues have been assigned (indicated by
     fit.nside[j]<0 = - rottype)   Then write standard pdb
First goes through fits and looks for fits that have both sequence
  and segment assigned:
     For  each segment and write output pdb in order of a.a. number 
  A TER statement is put at end of each segment.
Then writes fits with segments not assigned.
  A TER statemente is put at end of each fit.
Then write fits with sequence not assigned.
  A TER statemente is put at end of each fit.
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
    4+fit[i].nside[j] is side chain atom number
*/
{
   FILE *fp;
   int ifit,resnum,atomnum,iseg,writeter;
   char segchar[5],reschar[20];

printf("Converting  %d fits to PDB format and write to file %s\n",fitnum[0],fitpdbfile);
   fp= fopen(fitpdbfile,"w");
   for(ifit=1;ifit<=fitnum[0];ifit++){ // write fits sequence  and segment assignment
       if( fit[ifit].val2 < -900 ){ // segment assigned
          iseg =  -fit[ifit].val2 - 900;
          itoa(iseg,segchar);
          atomnum = 0;
          if(fit[ifit].val1<0)
                resnum = -fit[ifit].val1 -1;//Need to subtact 1 because writeonefit start at resnum+1
          else
                resnum = 0;
          if(ifit==fitnum[0]) 
             writeter = 1;  // indicates write TER at end of fit
          else if( (-fit[ifit+1].val2 -900) != iseg)
             writeter = 1;
          else
             writeter = 0;
          writeonefittopdb(fp,den1,geo,fit,fitnum,fitpdbfile,ifit,reschar,
                 &resnum,&atomnum,segchar,writeter);// return resnum = last residue written
       } // segment assigned
       else{
          printf("ERROR IN  writefittopdbassign() - no segment assigned for fit %d\n",ifit);
       }
   }//ifit
   fprintf(fp,"END\n");
   fclose(fp);  
}
#endif

void writefittopdbgen(struct griddata& den1,struct geometry& geo,
            struct pepfit fit[NFIT],int fitnum[NFIT],
            char fitpdbfile[])
/*  If first fit has segment assigned (assume all have segment assigned)
      - call writefittopdbassign() 
  Else - call writefittopdb() - and write each  fits as different segment. 
*/
{     
    if(fitnum[0]>=1){
       if( fit[1].val2 < -900 ){ // segment assigned
         writefittopdbassign(den1,geo,fit,fitnum,fitpdbfile);// This writes sequence assingnd fits in order
       }
       else{
          writefittopdb(den1,geo,fit,fitnum,fitpdbfile); // Use this version if segment not assigned
       }
    }
    else{
      printf("ERROR - NO FITS IN FITFILE\n");
    }
}
