// Initial values called from main
#include "tordyn.h" // routines for torsiondynamics
#include "peptide.h"  // routines for peptide creation and drawing
#include "drawwindow.h"  // also calls global.h
#include "trace.h"
#include "sphere.h"  // routines for sphere pull down
#include "math.h"
#include "bone.h"
#include "pdbtolev.h" // define NATOMS, needed for assignseq.h
#include "assignseq.h"
#include "initialize.h"
#include "string.h"
#include "stdlib.h" 

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

void initgeo(struct geometry& geo,struct maiddatafile& maidfile,struct rotomer rot[])
{
   strcpy(geo.setfile,"set2");  // Ranked setnumber that is being written and read
   setup(geo); // and next two calls set up geometry constants for creating  peptides
   setuprotomer(rot); 
   findconstants(geo);
   assignlevconst(geo);  // Needed for rotomer routines
   geo.movecb=0;  // if = 0 do not move Cb in Powell refinement of rotomer
   geo.ringang=0; //if = 0 only use the default ring angle, else try others
   geo.maxcbmove = 0.5;
   geo.maxfusedist = 1.5;
   geo.goodfusedist = 1.5; // NOT USED anymore
   geo.fitn = 1; // As default extend at end of fit1
   geo.useset = 1; // Means use set in tracetobr  (for checkifbr() do not want to use set to see fit meet another fit
   geo.setval = -1; // default val indicating point set
   geo.setside = -1; // set value for side atoms  
   geo.setmain = -1; // set value for mainchain
   geo.userank = 0; // 0 - use all fits in fitfile; 1 = previous ranked fits; 2= use the current picked geo fit
   geo.nmaxext = 7;  // Maxinum of  residues extended at one end
//   geo.maxenddist = 3.5;//Used up to 11feb99 maximum distance Ca ends allowed to move during "bestden" rfinemnt
   geo.maxenddist = 5.5;//USED UP TO 29NOV99: maximum distance Ca ends allowed to move during "bestden" rfinemnt
   geo.maxenddist = 2.0;//Used to 5jun00 TRY SMALLER VALUE maximum distance Ca ends allowed to move during "bestden" rfinemnt
   geo.maxenddist = 2.5;//Tried after 5jun00 TRY Larger VALUE maximum distance Ca ends allowed to move during "bestden" rfinemnt
   geo.phenylden= 0; // if = 1 then find density only over the last rigid ring of a.a. reisdue - e.ge. phenylaline
   geo.nowset = 0; // = 1 if set  array is allocated
   geo.frpoor = 0.395; // Allowed fraction of poor atoms in fit
   geo.oneside = 0; // default - run standard dynnaicms on all atoms (if = 5, then just dynamics on sidechain #5)
   geo.newfit =0; // default - start tracefit over, do not use old fits
//   geo.newfit =1; // Continue tracefit - using fits in fitfile
   geo.itype = -1;
//   geo.frden = 0.4;  // if aveden of last residue < frden of prev aveden, then quit
   geo.frden = 0.6;  // if aveden of last residue < frden of prev aveden, then quit
//   geo.arange = 2.0; // if last atom with arange (angstroms) of boundary, then stop adding residues
   geo.arange = 3.0; // if last atom with arange (angstroms) of boundary, then stop adding residues
   geo.arange = 4.0; // if last atom with arange (angstroms) of boundary, then stop adding residues
   // Ranges and step sizes for global (rigid) maximazation routine
   geo.thtot = PI/6.0; //range over which theta is varied
   geo.phitot = 2.0*PI; 
   geo.fixend = 0;  // Default - moves all atoms in tordyn - probably not need tobe set here becaue is set before dyn called
   geo.ncbang=20; // # of rotation steps to find best cb2
   geo.endmove = 4; //Value used most of 1998 - number of residues to move during end tor. dyn fit
//   geo.endmove = 6; //Default  - number of residues to move during end tor. dyn fit
//   geo.endmove = 60; //NOTE: this is now set in each call to initialfit or from menuThis will make tordyn use all residues while extending
#if 0 // works fairly well 
   geo.trandel = 0.75; // translatino step in angstroms for global minimizaton
   geo.tstep = 2; // number of translation steps  (max tran = tstep*trandel)
   geo.rotdel = 1.5; // maximu distance furthers point move during a rotation step
#endif
#if 1  // This keeps it within 0.6 A. of initial Ca position
   geo.trandel = 0.3; // translatino step in angstroms
   geo.tstep = 2; // number of translation steps  (max tran = tstep*trandel)
//   geo.tstep = 3; // number of translation steps  (max tran = tstep*trandel)
   geo.axstep = 1; // step + or -1 large step of geo.axdel along Ca to Ca axis
   geo.axdel = 1.2;  // = Angstrom step along axis
   geo.rotdel = 1.0; // maximu distance furthers point move during a rotation step
#endif
#if 1 // This change improves fit for at least one sheet (11feb99)
   geo.thtot = geo.thtot/2.0;
   geo.rotdel = geo.rotdel/2.0;
#endif
   geo.nstartres =3; // Number of residues in initial trial
   geo.nend[0]=geo.nend[1]=0; // set so 0 length
   geo.minrigid = 0.0; // If atom below this value - then added to nbad
   geo.poorrigid = 0.5*maidfile.denscale;//Used to count poor values
   geo.npoor = 1000; //Number of npoor atoms allowed - Default very large so not limiting
#if 1 // standard version - works fairly well
   geo.nbadend = 1; //Used before 1oct98  Max # of allowed bad (den < minrigid) for terminal residue
//   geo.nbadend = 0; // Max # of allowed bad (den < minrigid) for terminal residue
//   geo.usecb = 1;  //Used before March 5, 99: geo.usecb == 1, then use cb denisty  in avectermden and aventermden(Never use Cb2)
   geo.usecb = 2;  //Used after March 5, 99: geo.usecb == 1, then use cb denisty  in avectermden and aventermden(Never use Cb2)
#endif
   geo.nbadwhole = 0; // Max # of allowed bad (den < minrigid) for initial peptide
//   geo.nbadwhole = 1; // Max # of allowed bad (den < minrigid) for initial peptide
   geo.nbadwt = 0.5*maidfile.denscale; //subtract nbadwt*nbad/residue from avedensity 
   geo.savlen = 6;// USED up to 3dec99:After extension, total length must be this long for fit to be saved
//   geo.savlen = 5;//Try  shorter lenght -getsome bad fits for especially for sheets
   geo.extendden = 1.0*maidfile.denscale; //Used up to 30may00 ave density after global refinement must be this large before extended, also used to weight density error function
//   geo.extendden = 0.8*maidfile.denscale; //Tried on 30may00 ave density after global refinement must be this large before extended, also used to weight density error function
   geo.sideden = 0.7*geo.extendden; //Used most of time  Density for extendign side chain
//   geo.sideden = 0.5*geo.extendden; // Density for extendign side chain
   geo.angweight = 0.5; // This value works well   weight for angle error in density error function  
//   geo.angweight = 0.1; //This cause trouble    weight for angle error in density error function 
//   geo.angweight = 1.0; // weight for angle error in density error function 
   geo.bondwt = 1.0; // weight for error in Ca bond lenght
//   geo.setrad= 2; //OLD used before 16feb99 number of grid points about atom for set[][][] = -1
   geo.setrad= 1; //New used after 16feb99 number of grid points about atom for set[][][] = -1
   geo.spherelim=1; //if = 1 In trace.c++, only use start points that are in spheres in spherfile
   geo.nside = 2;  // default value
   geo.torend =1; //NOTE: version of avectermden and avaentermden that use this have been commented out Number of residues whose density is averaged in avectermden (=1 is default, just end residue)
//   geo.maxside = 4; // maximum number of side chains
   geo.maxside = 5; // maximum number of side chains
//   geo.minavecb2 = 0.5*maidfile.denscale;  // Minimum allowed value of ave Cb2 density to accempt final fit
   geo.minavecb2 = geo.minavecb= maidfile.denscale;  //Value used in avecb2 is Cb  density to accempt final fit
//   geo.minfrbad= 0.5;  // Used before Dec 8,98Minimum allowed fraction of residues allowed to have bad atoms
//   geo.minfrbad= 0.25;  // Used after Dec 8,98Minimum allowed fraction of residues allowed to have bad atoms
   geo.minfrbad= 0.3;  // Used after Dec 9,98  Need bigger values for sheets in shad
   geo.minfrbad= 0.2;  // Used after Dec 1,98  
   geo.curdist[1] = -1100.0; // used as indicator that trace has not been called and need to allocate set
   geo.ntypes = 7; //The default number of constaint types :Type 1 = helix, 2 = sheet, 3  = left helix 4= epsilon 5 = delta  6= gammar 7 = gammal 
//   geo.ca1_4min = 5.0; //Used upto 30may00: Minimum Ca[i] to Ca[i+4] dist for helix
   geo.ca1_4min = 4.7; //NEW - 30may00:  Minimum Ca[i] to Ca[i+4] dist for helix
   geo.linewidth = 1.0; 
   geo.fuse[0]= -1;  //fuse[0] = -1 (default); 0,1 or 2 (number of fits picked); fuse[1,2] = nterm, cterm fitnumber picked;  
   geo.fuse[1]=geo.fuse[2]= -1; // 1: Nterm, 2: cterm  if = -1, not picked yet
   geo.seqnum = -1; // indicates have not yet read sequence file
   geo.fnum = -1; // Number of geo residue that was picked
   geo.removefit = -1; // Do not remove picked fit
   geo.main[0] = 0; // Indicates first time called extend fit
   geo.checkmaindist = 0; // Indicates do not check distance from edited main in tordylnbat
   geo.lmax =2.0; // maximum distance from atom center that extend atoms density
   geo.tempfact = 60;
   geo.denmult = -100000.0;
#if 1 // New 19jun00 - now defined globally
    geo.valueweight[0] = 0.3; // value = -(aveden - maidfile.denscale*(valueweight[0]*nbad + valueweight[1]*badset) );
    geo.valueweight[1] = 1.0;
#endif
   geo.maxngap = 10;// Maximum allowed baseline (den0.8) length of edited gap connection to be used for extending
   geo.gapnum=1; // gap used in connectnum() varies from cnum-geo.gapnum to cnum+geo.gapnum
//   geo.rankmax = 25.0; // best ranking must be at least this good if gap OK
   geo.rankmax = 15.0; // New - after 4oct00 -best ranking must be at least this good if gap OK
   geo.rankgap = 5.0;//Difference between first and second best if use gap.
   geo.gnmax = 10; // (=tordata.nmax) default value - NOTE reset in extendfit.c++ based on resolution 
   geo.maxline = 1.0;  //Set default to very large value so only used for short fits-- maximum allowed values for short fits of ave line dist in bonedistgeo
   geo.minnres = 13; //Minimum number of resdues required to try and test gap in extendfit()
   geo.sidemaindist =2.5; // Minimum allowed distance for long sidechain atoms to mainchain 
   geo.miniden=0; //0 corresponds to boneden  = 0.8  -  the minimum value of bone iden used when looking for connection
   geo.addseed=0;// add to ipt when set srand48
} 

void initialize(struct geometry& geo,int fitnum[NFIT],
       struct pepfit fit[NFIT],struct spheredata sp[NSPHERES],
       struct maiddatafile& maidfile,struct tordyn& tordata,
       struct tordyn& tordata2,struct bonedata bone[NMAPS],
       struct pepset fitset,struct rotomer rot[])
{
   int i;
    maidfile.denscale=20.0; /*multiplier of floating point density before convert to
         type char.  In formatted map den=6 is 6 x stand. dev above aveage.
         6x20=120, anything above (maxchar-1)=126 set to 126; must keep negative density
         so gradient calculated correctly*/
    maidfile.nsphere = 0; // The total number of sphere
    maidfile.start = 1;  // This indicates just starting - draw simple cube
    fitnum[0]=0; // Indicates that there are no save geofits
    for(i=1;i<NFIT;i++){ //(This is repeated in tracepep) set fitnum < 0 - indicate is empty
      fitnum[i]= -10;
      fit[i].nalloc = -10;
   }
   for(i=0;i<NSPHERES;i++){ //Initially, do not move or pick sphere
      sp[i].picksphere=sp[i].movesphere= -1;
      sp[i].radius = 5.0; // initial radius in angstroms
   }
   setupauto(geo,tordata);
   setupauto(geo,tordata2);
   initgeo(geo,maidfile,rot);
   for(i=1;i< NMAPS;i++)
       bone[i].minden = 1.4;  //Default minimum value of density used in bone calculation
   fitset.rank[0][0]=0;
}

#if 1 // moved from menu.c++
void writemaidfile(int inew,struct maiddatafile& maidfile,struct griddata map[],
         struct bonedata bone[])
/* if inew = 1; then new - must create default file names,
    if = 0, then old, just use already created names
*/
{ 
    char ich[10],pdbname[10],denname[10],bonename[10];
    int i;
    FILE *fp,*tempfp; 

    printf("Writing maid file %s\n",maidfile.name);
    fp= fopen(maidfile.name,"w");
    fprintf(fp,"Number of pdbfiles:  %d\n",maidfile.npdb);
    for(i=1;i<=maidfile.npdb;i++){
          itoa(i,ich);
          strcpy(pdbname,"pdbfile");
          strcat(pdbname,ich); 
          fprintf(fp,"%s:  %s\n",pdbname,maidfile.pdb[i].name);
    }    
    fprintf(fp,"Number of denfiles:  %d\n",maidfile.nden);
    for(i=1;i<=maidfile.nden;i++){
         itoa(i,ich);
         strcpy(denname,"denfile");
         strcat(denname,ich);
         fprintf(fp,"%s:   %s\n",denname,map[i].mapname);
#if 0
         strcpy(corename,"corefile");
         strcat(corename,ich);
         core[i].num = i; // Indentity number of corefile
         fprintf(fp,"%s:   %s\n",corename,core[i].corename);
#endif
#if 1
         strcpy(bonename,"bonefile");
         strcat(bonename,ich);
         bone[i].num = i; // Indentity number of bonefile
         fprintf(fp,"%s:   %s\n",bonename,bone[i].bonename);
#endif
    }
    if(inew==1){ // new - create spherefile and write default value to it and create other files
       strcpy(maidfile.spherefile,maidfile.name);
       strcat(maidfile.spherefile,"spherefile");
       tempfp = fopen(maidfile.spherefile,"w");
       fprintf(tempfp,"%d\n",-1);// Input default  = -1 spheres - indicates do not use speres
       fclose(tempfp);
       strcpy(maidfile.fitfile,maidfile.name);
       strcat(maidfile.fitfile,"fit"); 
       strcpy(maidfile.pdbsequence,maidfile.name);
       strcat(maidfile.pdbsequence,"pdbsequence"); 
       strcpy(maidfile.editbones,maidfile.name);
       strcat(maidfile.editbones,"editbones"); 
       strcpy(maidfile.nsym,maidfile.name);
       strcat(maidfile.nsym,"symrottran"); 
       strcpy(maidfile.selenomet,maidfile.name);
       strcat(maidfile.selenomet,"selenomet"); 
    } 
    fprintf(fp,"%s:  %s\n","spherefile",maidfile.spherefile);
    fprintf(fp,"%s:  %s\n","fitfile",maidfile.fitfile);    
    fprintf(fp,"%s:  %s\n","editbonesfile",maidfile.editbones);    
    fprintf(fp,"%s:  %s\n","pdbsequencefile",maidfile.pdbsequence);    
    fprintf(fp,"%s:  %s\n","symrottranfile",maidfile.nsym);    
    fprintf(fp,"%s:  %s\n","selenometfile",maidfile.selenomet);    
    fclose(fp);
}

void readmap(char filename[],struct griddata& den1,
       struct maiddatafile& maidfile)
// New 12jul00 - must be sure can read long enough line
//    increase dimension to 2000
{
	char name[2000],order[3];
	int i,j,k,ntitle,nless,nmore;
	float dd,dds;
        FILE *fp;

        if( (fp=fopen(filename,"r")) == NULL){
          printf("error opening file %s\n",filename);
          exit(0);
        } 
        printf("Reading XPLOR map file %s\n",filename);
	fgets(name,200,fp); /* read top blank line */
	fgets(name,200,fp); /* read number of title lines */
	sscanf(name,"%d",&ntitle);
	for(i=1;i<=ntitle;i++)
             fgets(name,1999,fp);  /* read title - discard */
	fgets(name,200,fp); /* read section data */
printf("%s\n",name);
	sscanf(name,"%d%d%d%d%d%d%d%d%d",&den1.na[0],&den1.amin[0],
                       &den1.amax[0],&den1.na[1],&den1.amin[1],
                       &den1.amax[1], &den1.na[2],&den1.amin[2],
                       &den1.amax[2]);
printf("den1.na[0]=%d  den1.amin[0]=%d amax[0]=%d  na[1]=%d amin[1]=%d\n",
   den1.na[0],den1.amin[0],den1.amax[0],den1.na[1]); 
#if 0
        for(i=0;i<3;i++)
           den1.griddim[i]=den1.amax[i]-den1.amin[i]+2;/* the +2 is needed for iden array*/
        printf("Allocating den array dimensions %d %d %d\n",
                           den1.griddim[0],den1.griddim[1],den1.griddim[2]);
        scmat3(den1.den,den1.griddim[0],den1.griddim[1],den1.griddim[2]);// allocates array 3 dim array (mat)
                     // of type  unsigned char - n1,n2,n3,n4 size of each dimension
#endif
	fgets(name,200,fp); /* read unit cell data */
	sscanf(name,"%f%f%f%f%f%f",&den1.cell[1],&den1.cell[2],&den1.cell[3],
                                   &den1.cell[4],&den1.cell[5],&den1.cell[6]);
        initializemap(den1);
//        initializemapscreen(den1);
	fgets(name,200,fp); /* read matrix order */
	sscanf(name,"%3s",order);
	if( (order[0] != 'Z')  || (order[1]!='Y')  || (order[2]!='X') ){
            printf("ERROR: Not a formated XPLOR map\n");
            exit(0);
        }
	/* read matrix map values */
        nless=nmore=0;
	for(k=0;k<=den1.amax[2]-den1.amin[2];k++) {
                printf("Reading k section %d\n",k);
		fscanf(fp,"%s",name); /* read k section name and discard*/
		for(j=0;j<=den1.amax[1]-den1.amin[1];j++) {
	           for(i=0;i<=den1.amax[0]-den1.amin[0];i++) {
                      fscanf(fp,"%12f",&dd);
                      dds=dd*maidfile.denscale;
                      if(dds>MAXCHAR){//maxchar = 127 is maximum real density                         
                         dds=MAXCHAR;
                         nmore=nmore+1;
                      }
                      if(dds< -MAXCHAR+1){ // save dds = -MAXCHAR for use in bones routines
                           nless=nless+1;
                           dds= -MAXCHAR+1;
                      }
                      den1.den[i][j][k]=dds;
		     }
		}
	}
        printf("# of den > maxchar = %d; # < -maxchar = %d\n",nmore,nless);
	fclose(fp);
        printf("Done\n");
}
#endif

void setup(struct maiddatafile& maidfile,struct griddata map[],
          struct bonedata bone[],char name[])
/* Only required input is the name of the maid file
  Program outputs the maidfile name.maid which
    can be used as input for batch or graphics.
*/
{
  char mapname[300],pdbname[300],yesno[10],ich[10],ch[20];

   strcpy(maidfile.name,name);
   strcat(maidfile.name,".maid");
   printf("Input the name of map file\n");
   printf("    in NORMALIZED (mean=0,std=1) XPLOR or CNS text format:\n");
   scanf("%s",mapname);
   maidfile.nden=1;
   readmap(mapname,map[1],maidfile);// read the this map into struct griddata map[1]....
   strcpy(map[1].mapname,mapname);
   strcat(map[1].mapname,".den");
   itoa(1,ich);
   strcat(map[1].mapname,ich);
   saveden(map[1]);
   strcpy(bone[1].bonename,"none"); // indicates no bone file created yet
   printf("Is there a tentative PDB file (used only for comparison in graphics version)?\n");
   printf("   Note: Graphics MAID will center on the PDB- Enter yes/no\n");
   scanf("%s",yesno);
   if( (strcasecmp(yesno,"yes")==0)||(strcasecmp(yesno,"Y")==0) ){ //then read pdb file
       maidfile.npdb = 1;
       printf("Input name of PDB file in standard format\n");
       scanf("%s",pdbname);
       strcpy(maidfile.pdb[1].name,pdbname);
   }
   else if( (strcasecmp(yesno,"no")==0)||(strcasecmp(yesno,"N")==0) ){ //then read pdb file
      printf(" No pdb file input\n");
   }
   else{
      printf("ERROR: you must enter yes or no - quitting\n");
      exit(0);
   }
   printf("Input the skeletonization cut off density used for finding sheets and helices\n");
   printf("   For a good quality map use enter '1.4' - poorer try '1.2'\n");
   scanf("%f",&bone[1].minden);

   strcpy(bone[1].bonename,map[1].mapname);
   strcat(bone[1].bonename,".bone"); // bone  name is mapname with .bone_minden appended
   sprintf(ch,"%1.1f",bone[1].minden);
   strcat(bone[1].bonename,ch);
   findbone(map[1],bone[1]);
   writemaidfile(1,maidfile,map,bone);// "1" indicate create new file names
}    
   
#if 1 // moved from menu.c++
void saveden(struct griddata den1)
{
   long i,j,k;
   FILE *fps;

   printf("Saving the density in file %s\n",den1.mapname);
   fps=fopen(den1.mapname,"wb");
   fwrite(&den1.na[0],sizeof(den1.na[0]),3,fps);
   fwrite(&den1.amin[0],sizeof(den1.amin[0]),3,fps);
   fwrite(&den1.amax[0],sizeof(den1.amax[0]),3,fps);
   fwrite(&den1.cell[1],sizeof(den1.cell[1]),6,fps);
   for(k=0;k<=den1.amax[2]-den1.amin[2];k++) 
       for(j=0;j<=den1.amax[1]-den1.amin[1];j++) 
	    for(i=0;i<=den1.amax[0]-den1.amin[0];i++) 
                fwrite(&den1.den[i][j][k],sizeof(den1.den[i][j][k]),1,fps);
   fclose(fps);
   printf("Done\n");
}
#endif
