/* EGAD: minimization.cpp
   
Navin Pokala and Tracy Handel
Dept. of Molecular and Cell Biology
University of California, Berkeley
Copyright (C) 2003 Regents of the University of California
GNU Public License
Aug 12 2003

Absolutely no warranties are made or are implied with the use of this program or its parts.

This file contains functions for performing powell minimization of protein->Template in dihedral space.
Also contains functions for adjusting backbone torsions to compensate for distortions due to rebuilding 
	backbone structures with ideal geometry. 

The functions minimization_objective_function and convert_delta_dihedral_array_to_CHROMOSOME may also be
	used with other continuous dihedral optimization algorithms.

*/


#include "minimization.h"

int NUM_FUNCTION_CALLS;

/* given a base chr and delta_dihedral_array, generate chr_new (allocated by calling function).
	if  backbone_flag == 1, bkbnGENES are built.
	fixed_res!=NULL -> fix dihedrals for these seq positions. 
	float_res!=NULL -> float dihedrals only for these seq positions.
	delta_dihedral_array - all sidechain rotamer perturbations listed first, 
	then all backbone dihedral perturbations.
		ligands, pro, gly, cyd, ala sidechain dihedrals are defined as fixed,
		phi, psi, and omega are allowed to vary for backbone.
	chr_new can be used for CHROMOSOME_to_mini_pdbATOM 
*/

void convert_delta_dihedral_array_to_CHROMOSOME(CHROMOSOME *chr_new, double *delta_dihedral_array, CHROMOSOME *chr, 
											int backbone_flag, int *fixed_res, int *float_res)
{
  int dihed_ctr, q, fix_ctr, float_ctr, float_flag;
  
  dihed_ctr = 1;
    
  chr->genes = chr->firstgene;   
  chr_new->genes = chr_new->firstgene;
  fix_ctr=1; float_ctr=1; 
  while(chr->genes->nextgene!=NULL)
    {
	float_flag = 1;	/* by default, all residues can float */

	if(float_res != NULL)	/* if this is true, by default, only specified residues can float */
	{					
		float_flag = 0;	/* assume that this residue doesn't float, unless specified otherwise */
	
		while(float_res[float_ctr] < chr->genes->seq_position && float_res[float_ctr]!=ENDFLAG)
			++float_ctr;
	     	if(float_res[float_ctr] == chr->genes->seq_position)	
			float_flag = 1;

	}
	if(fixed_res != NULL)	
	{
	     while(fixed_res[fix_ctr] < chr->genes->seq_position && fixed_res[fix_ctr]!=ENDFLAG)
		++fix_ctr;
	     if(fixed_res[fix_ctr] == chr->genes->seq_position)	
			float_flag = 0;
	}
	
	/* keep CYD and ligands fixed; ALA is fixed by definition */
	if(chr->genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'c' || 
		chr->genes->choice_ptr->resparam_ptr->ligand_flag==1 ||
		chr->genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'A' ||
		chr->genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'G' ||
		chr->genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'P')
			float_flag=0;

      for(q=1;q<=chr->genes->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++q)
	{
  	  if(float_flag == 1)	/* move dihedrals only for residues that are allowed to float */
	  {
	  	chr_new->genes->chi[q] = chr->genes->chi[q] + delta_dihedral_array[dihed_ctr];
		++dihed_ctr;
	  }
	  else
		chr_new->genes->chi[q] = chr->genes->chi[q];
	  
	  
	  if(chr_new->genes->chi[q]>180)
	    chr_new->genes->chi[q] = chr_new->genes->chi[q] - 360.0;
	  else
	    if(chr_new->genes->chi[q]< -180)
		  chr_new->genes->chi[q] = 360 + chr_new->genes->chi[q];
	}
      chr->genes = chr->genes->nextgene;
      chr_new->genes = chr_new->genes->nextgene;
    }
   chr->genes = chr->firstgene;
   chr_new->genes = chr_new->firstgene;
   
  if(backbone_flag == 1)
  {
  chr->bkbngenes = chr->first_bkbngene;
  chr_new->bkbngenes = chr_new->first_bkbngene;
	fix_ctr=1; float_ctr=1; 
  while(chr->bkbngenes->nextbkbngene!=NULL)
    {
	float_flag = 1;	/* by default, all residues can float */

	if(float_res != NULL)	/* if this is true, by default, only specified residues can float */
	{					
		float_flag = 0;	/* assume that this residue doesn't float, unless specified otherwise */
		
		while(float_res[float_ctr] < chr->bkbngenes->seq_position && float_res[float_ctr]!=ENDFLAG)
			++float_ctr;
	     	if(float_res[float_ctr] == chr->bkbngenes->seq_position)	
			float_flag = 1;

	}
	if(fixed_res != NULL)	
	{
	     while(fixed_res[fix_ctr] < chr->bkbngenes->seq_position && fixed_res[fix_ctr]!=ENDFLAG)
		++fix_ctr;
	     if(fixed_res[fix_ctr] == chr->bkbngenes->seq_position)	
			float_flag = 0;
	}

	if(float_flag == 1)	/* move dihedrals only for residues that are allowed to float */
	{
	  chr_new->bkbngenes->phi = chr->bkbngenes->phi + delta_dihedral_array[dihed_ctr];
	  ++dihed_ctr;
  	}
	else
	  chr_new->bkbngenes->phi = chr->bkbngenes->phi;

          if(chr_new->bkbngenes->phi>180)
	    chr_new->bkbngenes->phi = chr_new->bkbngenes->phi -360.0;
	  else
	    if(chr_new->bkbngenes->phi < -180)
		  chr_new->bkbngenes->phi = chr_new->bkbngenes->phi + 360.0;
				 

	if(float_flag == 1)
	{
		chr_new->bkbngenes->psi = chr->bkbngenes->psi + delta_dihedral_array[dihed_ctr];
		++dihed_ctr;
	}
    	else
		chr_new->bkbngenes->psi = chr->bkbngenes->psi;
	 
	 
	 if(chr_new->bkbngenes->psi>180)
	    chr_new->bkbngenes->psi = chr_new->bkbngenes->psi -360.0;
	  else
	    if(chr_new->bkbngenes->psi < -180)
		  chr_new->bkbngenes->psi = chr_new->bkbngenes->psi + 360.0; 


	if(float_flag == 1)
	{
		chr_new->bkbngenes->omega = chr->bkbngenes->omega + delta_dihedral_array[dihed_ctr];
		++dihed_ctr;
	}
	else
		chr_new->bkbngenes->omega = chr->bkbngenes->omega;
	
	 
	 if(chr_new->bkbngenes->omega>180)
	    chr_new->bkbngenes->omega = chr_new->bkbngenes->omega -360.0;
	  else
	    if(chr_new->bkbngenes->omega < -180)
		  chr_new->bkbngenes->omega = chr_new->bkbngenes->omega + 360.0; 

		  
	chr->bkbngenes = chr->bkbngenes->nextbkbngene;
	chr_new->bkbngenes = chr_new->bkbngenes->nextbkbngene;
    }
   chr->bkbngenes = chr->first_bkbngene;
   chr_new->bkbngenes = chr_new->first_bkbngene;
   }
}


/*  returns value of minimization_objective_function.
    delta_dihedral_array - defined as described for convert_delta_dihedral_array_to_CHROMOSOME.
    builds the structure and evaluates the energy E.
    if backbone is allowed to move (protein->parameters.rebuild_backbone_flag==1), and if RESTRAIN_MINIMIZATION_FLAG=1,
	the objective function is: if(E>=0), Fobj = E*rmsd; else Fobj = E/rmsd; 
		(rmsd = current structure-starting structure bkbn rmsd for N, CA, C)
    if RESTRAIN_MINIMIZATION_FLAG=0, E is returned.
*/

double minimization_objective_function(double *delta_dihedral_array, POWELL *powell)
{
    static mini_pdbATOM *pdb=NULL;
    static mini_pdbATOM *sideatoms=NULL;
    static double best_so_far;
    static int max_atoms=0;
    extern int TORSION_FLAG, COULOMB_FLAG, GB_FLAG, SASA_FLAG;
    double E_total, rmsd, r2;
    static char *filename;
    FILE *file_ptr;
    static SIDECHAIN side;
    var_fix_ENERGY varfix_energy;
    int i, r, q, j, float_res_ctr, include_flag, exclude_res_ctr, n;
    PROTEIN *protein;
    extern int LOGFILE_FLAG, RESTRAIN_MINIMIZATION_FLAG;
    extern double VDW_CUTOFF;

    protein = powell->protein;

    rmsd = 1.0; float_res_ctr=1;
    

	if(protein->parameters.max_iterations == ENDFLAG || protein->parameters.max_iterations ==2)
		return(0);	

    if(pdb == NULL)
	{
	    pdb = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
	    sideatoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
	    side.atom = (micro_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(micro_pdbATOM));
	    filename = (char *)calloc(MAXLINE,sizeof(char));
	    best_so_far = DBL_MAX;
	    max_atoms=MAX_ATOMS;
	}
    if(max_atoms<MAX_ATOMS)	/* MAX_ATOMS has gotten bigger; resize static arrays */
    {
	free_memory(pdb); free_memory(sideatoms);
	pdb = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
	sideatoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
	max_atoms=MAX_ATOMS;
    }
    
    convert_delta_dihedral_array_to_CHROMOSOME(&protein->chr[2], delta_dihedral_array, &protein->chr[1], protein->parameters.rebuild_backbone_flag, protein->min_fixed_res, protein->min_float_res);
    
	/* the backbone is allowed to move */
    if(protein->parameters.rebuild_backbone_flag == 1)
    {
	/* generate the structure and score the energy */
	CHROMOSOME_to_mini_pdbATOM(&protein->chr[2], protein->mini_fixed_atoms, pdb, protein->chain_anchor_bkbn);
	protein->chr[2].energy_list = energy_calc(pdb, GB_FLAG, SASA_FLAG);
	protein->chr[2].energy_list.E_1_4 = ((double)TORSION_FLAG)*protein->chr[2].energy_list.E_1_4;
	protein->chr[2].energy_list.E_coulomb = ((double)COULOMB_FLAG)*protein->chr[2].energy_list.E_coulomb;
    
		E_total = ENERGY_TOTAL(protein->chr[2].energy_list);


	if(RESTRAIN_MINIMIZATION_FLAG != 0)	/* the minimization is restrained to keep the structure close to the original structure */
    	{					
			
		if(RESTRAIN_MINIMIZATION_FLAG==1)
		{
			rmsd=0;	 float_res_ctr=1; exclude_res_ctr=1;
			protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
			protein->chr[2].bkbngenes = protein->chr[2].first_bkbngene;
			while(protein->chr[2].bkbngenes->seq_position!=ENDFLAG)
			{
	
				include_flag = 1;	// include all positions for restraint by default 
			
			if(protein->torsion_exclude_res!=NULL)	// these positions have been defined to be ignored for restraint 
			{
				while(protein->torsion_exclude_res[exclude_res_ctr] < protein->chr[2].bkbngenes->seq_position && 
						protein->torsion_exclude_res[exclude_res_ctr] != ENDFLAG)
					++exclude_res_ctr;
				if(protein->torsion_exclude_res[exclude_res_ctr] == protein->chr[2].bkbngenes->seq_position)
					include_flag = 0;
			}
			
			// square-well; could easily be replaced by a harmonic potential
			if(include_flag == 1)
			{
			if(fabs(protein->chr[2].bkbngenes->phi - protein->chr[1].bkbngenes->phi) > protein->parameters.max_bkbn_delta_dihedral)
					++rmsd;
			if(fabs(protein->chr[2].bkbngenes->psi - protein->chr[1].bkbngenes->psi) > protein->parameters.max_bkbn_delta_dihedral)
					++rmsd;
			if(fabs(protein->chr[2].bkbngenes->omega - protein->chr[1].bkbngenes->omega) > protein->parameters.max_bkbn_delta_dihedral)
					++rmsd;
			}
				protein->chr[1].bkbngenes = protein->chr[1].bkbngenes->nextbkbngene;
				protein->chr[2].bkbngenes = protein->chr[2].bkbngenes->nextbkbngene;

			}
			protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
			protein->chr[2].bkbngenes = protein->chr[2].first_bkbngene;

			E_total += VDW_CUTOFF*rmsd; 

		}
		
		if(RESTRAIN_MINIMIZATION_FLAG==2)  // score = (E>0) E_total*rmsd or (E<0) E_total/rmsd
		{
		i=1; j=1; float_res_ctr=1; exclude_res_ctr=1; n=0; rmsd = 0;		

 		while(protein->bkbn[i].seq_position!=ENDFLAG)
		{
			include_flag = 1;	/* include all positions for rmsd calc by default */

			
			if(protein->torsion_exclude_res!=NULL)	/* these positions have been defined to be ignored for rmsd */
			{
				while(protein->torsion_exclude_res[exclude_res_ctr] < protein->bkbn[i].seq_position && 
						protein->torsion_exclude_res[exclude_res_ctr] != ENDFLAG)
					++exclude_res_ctr;
				if(protein->torsion_exclude_res[exclude_res_ctr] == protein->bkbn[i].seq_position)
					include_flag = 0;
			}
			
			if(include_flag == 1)
			{
				while(pdb[j].seq_position!=protein->bkbn[i].seq_position)
					++j;
				
				r2 = Distance_sqrd(pdb[j].coord, protein->bkbn[i].N);  rmsd += r2; ++n; 	 ++j;
			/*	r2 = Distance_sqrd(pdb[j].coord, protein->bkbn[i].H);  rmsd += r2; ++n; */ ++j; 
				r2 = Distance_sqrd(pdb[j].coord, protein->bkbn[i].CA); rmsd += r2; ++n; 	 ++j; 
			/*	r2 = Distance_sqrd(pdb[j].coord, protein->bkbn[i].HA); rmsd += r2; ++n; */ ++j;
				r2 = Distance_sqrd(pdb[j].coord, protein->bkbn[i].C);  rmsd += r2; ++n; 	 ++j;
			/*	r2 = Distance_sqrd(pdb[j].coord, protein->bkbn[i].O);  rmsd += r2; ++n; */ ++j;
			/*	r2 = Distance_sqrd(pdb[j].coord, protein->bkbn[i].CB); rmsd += r2; ++n; */ ++j;
				
				while(pdb[j].seq_position == protein->bkbn[i].seq_position)
					++j;
			}
			++i;
		}
		rmsd = sqrt(rmsd/n);

	
		if(E_total>=0)
    			E_total = E_total*rmsd;
		else
		{
			if(rmsd>TINY)
				E_total = E_total/rmsd;
			else						/* prevent divide by 0 */
				E_total = E_total/TINY;
		}


		}
    	}	
    }
    else	/* only sidechains are moving */
	{
		/* Instead of using CHROMOSOME_to_mini_pdbATOM to generate the entire structure, 
				we can save a bit of time by dealing with two mini_pdbATOM arrays...keep reading */

		/* create sideatoms */
	    protein->chr[2].genes = protein->chr[2].firstgene;  
	    i=1;
	    while(protein->chr[2].genes->seq_position!=ENDFLAG)
	    {	
		sideatoms[i].coord = protein->chr[2].genes->choice_ptr->bkbn->CB;
		sideatoms[i].seq_position = protein->chr[2].genes->seq_position;
		sideatoms[i].sasa=0;
		sideatoms[i].atom_ptr = &protein->chr[2].genes->choice_ptr->resparam_ptr->atom[8];
		++i;
		
		build_a_SideChain(&side, protein->chr[2].genes->choice_ptr->bkbn,
				    protein->chr[2].genes->choice_ptr->resparam_ptr,&protein->chr[2].genes->seq_position, 
				    protein->chr[2].genes->chi); 

		r=9;
		for(q=1;q<=protein->chr[2].genes->choice_ptr->resparam_ptr->numberofSideAtoms;++q)
		{
		    sideatoms[i].coord = side.atom[r].coord; 

		    sideatoms[i].atom_ptr = &protein->chr[2].genes->choice_ptr->resparam_ptr->atom[r];
		    sideatoms[i].seq_position = protein->chr[2].genes->seq_position;
		    sideatoms[i].born_radius=0;
		    sideatoms[i].sasa = 0;
		    ++i; ++r;   
		}
		
		protein->chr[2].genes = protein->chr[2].genes->nextgene;
	    }
	    sideatoms[i].seq_position=ENDFLAG;
	    sideatoms[i].sasa=ENDFLAG;
	    sideatoms[i].atom_ptr=NULL;
	    protein->chr[2].genes = protein->chr[2].firstgene;  
	    
		/* internal energy of sidechain atoms and inter-sidechain energies */
	    protein->chr[2].energy_list = energy_calc(sideatoms, GB_FLAG, SASA_FLAG);

		/* energy between sidechain atoms and the fixed backbone atoms */
	    varfix_energy = var_fixed_energy_calc(protein->CBminus_minipdb, sideatoms,1);
	    
	/* this saves us from calculating the internal <fixed> backbone energies each time;
		the BORN and SASA energies will be off, but since minimization usually has these off,
		it shouldn't be too bad..... */

	    protein->chr[2].energy_list.E_vdw += varfix_energy.E_vdw;
	    protein->chr[2].energy_list.E_coulomb += ((double)COULOMB_FLAG)*varfix_energy.E_coulomb;
	    protein->chr[2].energy_list.E_1_4 += ((double)TORSION_FLAG)*varfix_energy.E_1_4;
	    protein->chr[2].energy_list.E_born += ((double)GB_FLAG)*varfix_energy.E_born;
	    protein->chr[2].energy_list.E_pol += ((double)GB_FLAG)*varfix_energy.E_pol;
	    protein->chr[2].energy_list.E_sasa = ((double)SASA_FLAG)*varfix_energy.E_sasa;

		E_total = ENERGY_TOTAL(protein->chr[2].energy_list);

	}	
    
  
    
    
    
    ++NUM_FUNCTION_CALLS;
    if(E_total < best_so_far)
	{
	    best_so_far = E_total;
		 
            if(LOGFILE_FLAG != 0)
	    if(NUM_FUNCTION_CALLS%50 == 0 || NUM_FUNCTION_CALLS==1)
            {
	    sprintf(filename, "%s.log", protein->parameters.output_prefix);
	    file_ptr = fopen_file(filename, "a");
	    fprintf(file_ptr, "energy eval %d\t%f\n", NUM_FUNCTION_CALLS, best_so_far);
	    fclose(file_ptr);
	    }
	}

    return(E_total);
}


/*
	For bkbn, an array of CARTESIAN pointers is assigned such that
		cartesian[i] = &bkbn.X (X=N,H,CA,etc).
	if exclude_res!=NULL, the listed positions are not included.
	if include_res!=NULL, only those listed positions are included.
	called by bkbn_torsion_overlay_rmsd.
*/
void BACKBONE_to_CARTESIAN_array(BACKBONE *bkbn, CARTESIAN **cartesian, int *exclude_res, int *include_res)
{
	int i,n;
	int include_flag, exclude_ctr, include_ctr;

	n=1; i=1; exclude_ctr = 1; include_ctr = 1;
	while(bkbn[n].seq_position!=ENDFLAG)	
	{	
		include_flag = 1; /* include by default */
		if(include_res != NULL)	/* if this is true, not included unless explicitly stated */
		{
			include_flag = 0;
			while(include_res[include_ctr] < bkbn[n].seq_position && include_res[include_ctr] != ENDFLAG)
				++include_ctr;
			if(include_res[include_ctr] == bkbn[n].seq_position)
				include_flag = 1;
		}
		if(exclude_res != NULL)	
		{
			while(exclude_res[exclude_ctr] < bkbn[n].seq_position && exclude_res[exclude_ctr] != ENDFLAG)
			{
				++exclude_ctr;
			}
			if(exclude_res[exclude_ctr] == bkbn[n].seq_position)
			{
				include_flag = 0;
			}
		}
				
		if(include_flag == 1)
		{						
		cartesian[i] = &bkbn[n].N;		
		++i;
		cartesian[i] = &bkbn[n].CA;
		++i;
		cartesian[i] = &bkbn[n].C;
		++i;
	 	cartesian[i] = &bkbn[n].O;
		++i;
		cartesian[i] = &bkbn[n].CB;
		++i;
		cartesian[i] = &bkbn[n].H;
		++i;
		cartesian[i] = &bkbn[n].HA;  
		++i;   
		}

		++n;
	}
	cartesian[i] = NULL;
}


/*
	for delta_dihedral_array, this function returns the rmsd (rmsd = current structure-starting structure bkbn rmsd).
	working_protein->torsion_exclude_res!=NULL --> listed positions not included in scoring.
	working_protein->torsion_include_res!=NULL --> only listed positions scored.
	delta_dihedral_array = array of backbone dihedral perturbations (phi,psi,omega) for all residues
*/
double bkbn_torsion_overlay_rmsd(double *delta_dihedral_array, POWELL *powell)
{	
	int dihed_ctr,  i,  n, first_time_flag, chain_number;
  	CHROMOSOME *chr;
	static BACKBONE *bkbn=NULL;
	static CARTESIAN **cartesian_star, **original_cartesian_star;
	static int max_residues=0;
	double dx,dy,dz,rmsd;
	div_t chain_id;
	PROTEIN *working_protein;

	working_protein = powell->protein;

	/*
	for conservation of inter-atomic distances

	static float *original_distance_list;
	extern int LOGFILE_FLAG;
	static char *filename=NULL;
	FILE *file_ptr;
	int j,q;
	double r;

	if(filename==NULL)
		filename = (char *)calloc(MAXLINE,sizeof(char));
	*/

        
	/* allocate memory */
	first_time_flag = 0;
	if(bkbn == NULL)
	{
		first_time_flag = 1;
		bkbn = (BACKBONE *)calloc(MAX_RESIDUES, sizeof(BACKBONE));
		cartesian_star = (CARTESIAN **)calloc(7*MAX_RESIDUES,sizeof(CARTESIAN *));
		original_cartesian_star = (CARTESIAN **)calloc(7*MAX_RESIDUES,sizeof(CARTESIAN *));
		max_residues=MAX_RESIDUES;
	}
	if(max_residues<MAX_RESIDUES) /* MAX_RESIDUES has gotten bigger; resize static arrays */
	{
		free_memory(bkbn); free_memory(cartesian_star); free_memory(original_cartesian_star);
		
		bkbn = (BACKBONE *)calloc(MAX_RESIDUES, sizeof(BACKBONE));
		cartesian_star = (CARTESIAN **)calloc(7*MAX_RESIDUES,sizeof(CARTESIAN *));
		original_cartesian_star = (CARTESIAN **)calloc(7*MAX_RESIDUES,sizeof(CARTESIAN *));
		max_residues=MAX_RESIDUES;
	}

	chr = &working_protein->chr[1];

	/* build backbone using base dihedrals in working_protein->chr[1] and 
		dihedral angle perturbations in delta_dihedral_array */
  	dihed_ctr = 1;
  chr->bkbngenes = chr->first_bkbngene;
  i=1; chain_number = 0; 
  while(chr->bkbngenes->nextbkbngene!=NULL)
    {
	  bkbn[i].seq_position = chr->bkbngenes->seq_position;

	  chain_id = div(bkbn[i].seq_position, 1000);

	  if(chain_number != chain_id.quot + 1)
	  {
		chain_number = chain_id.quot + 1;
		bkbn[i] = working_protein->chain_anchor_bkbn[chain_number];
	  }

	  bkbn[i].phi = chr->bkbngenes->phi + delta_dihedral_array[dihed_ctr];
	  ++dihed_ctr;
          if(bkbn[i].phi>180)
	    bkbn[i].phi = bkbn[i].phi -360.0;
	  else
	    if(bkbn[i].phi < -180)
		  bkbn[i].phi = bkbn[i].phi + 360.0;
				     
	 bkbn[i].psi = chr->bkbngenes->psi + delta_dihedral_array[dihed_ctr];
	 ++dihed_ctr;
	 if(bkbn[i].psi>180)
	    bkbn[i].psi = bkbn[i].psi -360.0;
	  else
	    if(bkbn[i].psi < -180)
		  bkbn[i].psi = bkbn[i].psi + 360.0; 

	
	   bkbn[i].omega = chr->bkbngenes->omega + delta_dihedral_array[dihed_ctr];
		++dihed_ctr;
          if(bkbn[i].omega>180)
	    bkbn[i].omega = bkbn[i].omega -360.0;
	  else
	    if(bkbn[i].omega < -180)
		  bkbn[i].omega = bkbn[i].omega + 360.0;

	    ++i;

		  
	chr->bkbngenes = chr->bkbngenes->nextbkbngene;
    }
   chr->bkbngenes = chr->first_bkbngene;
  
	bkbn[i].seq_position = ENDFLAG;	
	
	build_BACKBONE(bkbn);

/* 
	this commented-out region is for using conservation of inter-atomic distances as the objective function, 
	rather than the more conventional rmsd with the known input structure  
	 
	if(first_time_flag == 1)
	{
		original_distance_list = (float *)calloc(25*MAX_RESIDUES*(MAX_RESIDUES-1),sizeof(float));
	BACKBONE_to_CARTESIAN_array(working_protein->bkbn, original_cartesian_star, 
								working_protein->torsion_exclude_res, working_protein->torsion_include_res);
		n=1;
		i=1;
		while(original_cartesian_star[i]!=NULL)
		{
			if(original_cartesian_star[i]->x !=0)
			if(original_cartesian_star[i]->y !=0)
			if(original_cartesian_star[i]->z !=0)
			{
			j=i+1;
			while(original_cartesian_star[j]!=NULL)
			{
				dx = original_cartesian_star[i]->x - original_cartesian_star[j]->x; 
				dy = original_cartesian_star[i]->y - original_cartesian_star[j]->y; 
				dz = original_cartesian_star[i]->z - original_cartesian_star[j]->z; 
				original_distance_list[n] = dx*dx  + dy*dy + dz*dz;
				++j; ++n;
			}
			}
			++i;
		}
	}
	BACKBONE_to_CARTESIAN_array(bkbn, cartesian_star, working_protein->torsion_exclude_res, working_protein->torsion_include_res);
	rmsd = 0;
	n=1;
	i=1;
	while(cartesian_star[i]!=NULL)
		{
			if(original_cartesian_star[i]->x !=0)
			if(original_cartesian_star[i]->y !=0)
			if(original_cartesian_star[i]->z !=0)
			{
			j=i+1;	
			while(cartesian_star[j]!=NULL)
			{
				dx = cartesian_star[i]->x - cartesian_star[j]->x; 
				dy = cartesian_star[i]->y - cartesian_star[j]->y; 
				dz = cartesian_star[i]->z - cartesian_star[j]->z; 
				r = dx*dx  + dy*dy + dz*dz - original_distance_list[n];
				rmsd += r*r;
				++j; ++n;
			}
			}
			++i;
		}
	rmsd = sqrt(rmsd/n);
	*/

	
	/* traditional bkbn-bkbn rmsd between the original and ideal-geometry structure */	 

	if(first_time_flag == 1) /* calculate for the original protein once */
		BACKBONE_to_CARTESIAN_array(working_protein->bkbn, original_cartesian_star,
									working_protein->torsion_exclude_res, working_protein->torsion_include_res);
	BACKBONE_to_CARTESIAN_array(bkbn, cartesian_star, working_protein->torsion_exclude_res, working_protein->torsion_include_res);
	rmsd = 0;
	i=1; n=0;
	while(cartesian_star[i] != NULL)
	{
		if(original_cartesian_star[i]->x !=0)
		if(original_cartesian_star[i]->y !=0)
		if(original_cartesian_star[i]->z !=0)
		{
		++n;
		dx = cartesian_star[i]->x - original_cartesian_star[i]->x; 
		dy = cartesian_star[i]->y - original_cartesian_star[i]->y; 
		dz = cartesian_star[i]->z - original_cartesian_star[i]->z; 
		rmsd += dx*dx + dy*dy +dz*dz;
		}
		++i;
	}
	rmsd = sqrt(rmsd/n);
	 
	++NUM_FUNCTION_CALLS;
    	
	return(rmsd);
}


/* Manages JOBTYPE MINIMIZE and JOBTYPE IDEAL_BACKBONE_GEOMETRY
   For MINIMIZE, powell-minimize protein->Template in dihedral space 
	(backbone moves if protein->parameters.rebuild_backbone_flag==1). 
   For IDEAL_BACKBONE_GEOMETRY, rebuild backbone w/ ideal geometry, adjust torsions to minimize rmsd with 
	original structure.	
   Assumes protein has been initialized with input_stuff.
   final structure stored in protein->final_pdb for output_PROTEIN 
*/
void torsion_minimize_PROTEIN(PROTEIN *protein)
{
  int num_dihedrals;
  int i,j, stopflag;
  int float_ctr,fix_ctr, float_flag;
  VARIABLE_POSITION *dummy_varpos;
  pdbATOM *fixed_atoms;
  double **delta_dihedral_array=NULL, *energies=NULL;
  mini_pdbATOM  *mini_fixed_atoms;
  char *filename;
  FILE *file_ptr;
  extern int SASA_FLAG, LOGFILE_FLAG;
  extern CARTESIAN_VECTOR ***VECTOR_PTR;
  extern double MAX_OPTIMIZATION_TIME;
  time_t now, start_time;
  POWELL powell;
  

  powell.minimize_rotamer_datatype = NULL;
  powell.protein = protein;
  
  protein->sizeof_chr_array = 2;
  protein->chr = (CHROMOSOME *)calloc(3, sizeof(CHROMOSOME));
  
  if(SASA_FLAG == 1 && MAX_ATOMS <= REAL_MAX_ATOMS)
    {
      if(VECTOR_PTR==NULL)
	{   
	  VECTOR_PTR = (CARTESIAN_VECTOR ***)calloc(MAX_ATOMS, sizeof(CARTESIAN_VECTOR **));
	  for(i=0;i<MAX_ATOMS;++i)
	    VECTOR_PTR[i] = (CARTESIAN_VECTOR **)calloc(MAX_ATOMS, sizeof(CARTESIAN_VECTOR *));
	}
    }
  
  
  protein->CBminus_minipdb = NULL;
  
  if(LOGFILE_FLAG != 0)
    {
      filename = (char *)calloc(MAXLINE, sizeof(char));
      sprintf(filename, "%s.log", protein->parameters.output_prefix);
      file_ptr = fopen_file(filename, "w");
      fclose(file_ptr);
      free_memory(filename);
    }
  
  
  
  /* chr[1] is the Template; chr[2] is the working chromosome that will contain the modified dihedrals during minimization */    
  for(i=1;i<=2;++i)	
    pdbATOM_to_CHROMOSOME(protein->Template, &protein->chr[i], protein->resparam);
  
  
  mini_fixed_atoms = protein->mini_fixed_atoms;
  fixed_atoms = protein->fixed_atoms;
  
  protein->mini_fixed_atoms = NULL;
  protein->fixed_atoms = NULL;
  dummy_varpos = NULL;
  if(protein->parameters.rebuild_backbone_flag == 0)
    {
      dummy_varpos = (VARIABLE_POSITION *)calloc(MAX_RESIDUES, sizeof(VARIABLE_POSITION));
      protein->mini_fixed_atoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
      protein->fixed_atoms = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM));
      i=1; j=1;
      while(protein->Template[i].seq_position!=ENDFLAG)
	{
	  dummy_varpos[j].seq_position = protein->Template[i].seq_position;
	  while(dummy_varpos[j].seq_position == protein->Template[i].seq_position && protein->Template[i].seq_position!=ENDFLAG)
	    ++i;
	  ++j;        						    
	}
      dummy_varpos[j].seq_position = ENDFLAG;
      make_fixed_pdbATOM(protein->Template, dummy_varpos, protein->fixed_atoms);
      pdbATOM_to_mini_pdbATOM(protein->fixed_atoms,protein->mini_fixed_atoms); 
      free_memory(dummy_varpos);
      
      protein->CBminus_minipdb = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
      
      i=1; j=1;
      while(protein->mini_fixed_atoms[i].seq_position!=ENDFLAG)
	{
	  if(protein->mini_fixed_atoms[i].atom_ptr->atomnumber !=8)
	    {
	      protein->CBminus_minipdb[j] = protein->mini_fixed_atoms[i];
	      ++j;
	    }
	  ++i;
	}
      protein->CBminus_minipdb[j].seq_position=ENDFLAG;
      protein->CBminus_minipdb[j].sasa=ENDFLAG;
      protein->CBminus_minipdb[j].atom_ptr=NULL;
    }
  
  num_dihedrals = 0;
  if(strcmp(protein->parameters.algorithm, "ADJUST_BACKBONE_TORSIONS") != 0)
    {
      
      protein->chr[1].genes = protein->chr[1].firstgene;   
	float_ctr=1; fix_ctr=1;
      while(protein->chr[1].genes->nextgene!=NULL)
	{
										
	float_flag = 1;	/* by default, all residues can float */

	if(protein->min_float_res != NULL)	/* if this is true, by default, only specified residues can float */
	{					
		float_flag = 0;	/* assume that this residue doesn't float, unless specified otherwise */
	
		while(protein->min_float_res[float_ctr] < protein->chr[1].genes->seq_position && protein->min_float_res[float_ctr]!=ENDFLAG)
			++float_ctr;
	     	if(protein->min_float_res[float_ctr] == protein->chr[1].genes->seq_position)	
			float_flag = 1;

	}
	if(protein->min_fixed_res != NULL)	
	{
	     while(protein->min_fixed_res[fix_ctr] < protein->chr[1].genes->seq_position && protein->min_fixed_res[fix_ctr]!=ENDFLAG)
		++fix_ctr;
	     if(protein->min_fixed_res[fix_ctr] == protein->chr[1].genes->seq_position)	
			float_flag = 0;
	}
	
	/* keep CYD and ligands fixed; ALA is fixed by definition */
	if(protein->chr[1].genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'c' || 
		protein->chr[1].genes->choice_ptr->resparam_ptr->ligand_flag==1 ||
		protein->chr[1].genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'A' ||
		protein->chr[1].genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'G' ||
		protein->chr[1].genes->choice_ptr->resparam_ptr->one_letter_code[0] == 'P')
			float_flag=0;


	  if(float_flag==1)
	  	num_dihedrals += protein->chr[1].genes->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;


	  protein->chr[1].genes = protein->chr[1].genes->nextgene;
	}
      protein->chr[1].genes = protein->chr[1].firstgene;
      
      if(protein->parameters.rebuild_backbone_flag==1)
	{   
	  fprintf(stderr,"WARNING The backbone will be rebuilt using ideal bond geometry.\n");
	  fprintf(stderr,"\tThe initial template pdb file should have ideal bond geometry.\n");
	  fprintf(stderr,"\tIf it does not, rebuild it with JOBTYPE IDEAL_BACKBONE_GEOMETRY\n");	
	  fprintf(stderr,"\tOtherwise, the structure will be grossly perturbed.\n");

	  float_ctr=1; fix_ctr=1;

	  protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
	  while(protein->chr[1].bkbngenes->nextbkbngene!=NULL)
	    {

		float_flag = 1;	/* by default, all residues can float */

		if(protein->min_float_res != NULL)	/* if this is true, by default, only specified residues can float */
		{					
			float_flag = 0;	/* assume that this residue doesn't float, unless specified otherwise */
	
			while(protein->min_float_res[float_ctr] < protein->chr[1].bkbngenes->seq_position && 
					protein->min_float_res[float_ctr]!=ENDFLAG)
					++float_ctr;
	     		if(protein->min_float_res[float_ctr] == protein->chr[1].bkbngenes->seq_position)	
				float_flag = 1;
		}
		if(protein->min_fixed_res != NULL)	
		{
	     		while(protein->min_fixed_res[fix_ctr] < protein->chr[1].bkbngenes->seq_position && 
					protein->min_fixed_res[fix_ctr]!=ENDFLAG)
				++fix_ctr;
	     		if(protein->min_fixed_res[fix_ctr] == protein->chr[1].bkbngenes->seq_position)	
				float_flag = 0;
		}

	  	if(float_flag==1)
	  		num_dihedrals += 3;

	      	protein->chr[1].bkbngenes = protein->chr[1].bkbngenes->nextbkbngene;
	    }
	  protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
	} 
    }
  else
    {
      protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
      while(protein->chr[1].bkbngenes->nextbkbngene!=NULL)
	{
	  num_dihedrals += 3;
	  protein->chr[1].bkbngenes = protein->chr[1].bkbngenes->nextbkbngene;
	}
      protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
      
    }	    
  
  
  if(protein->parameters.max_iterations == 0)
    protein->parameters.max_iterations = num_dihedrals;
  if(protein->parameters.max_function_calls == 0)
    protein->parameters.max_function_calls = num_dihedrals*num_dihedrals*protein->parameters.max_iterations;
  
  
  
  energies = (double *)calloc(num_dihedrals+3, sizeof(double));


  delta_dihedral_array = (double **)calloc((num_dihedrals+3), sizeof(double *));
  for(i=0;i<(num_dihedrals+3);++i)
    delta_dihedral_array[i] = (double *)calloc((num_dihedrals+3), sizeof(double));
  
  NUM_FUNCTION_CALLS=0;
  
  for(j=1;j<=num_dihedrals;++j)
    delta_dihedral_array[1][j] = 0.0;
  energies[1] = 0;
  
  if(protein->final_pdb == NULL)
    protein->final_pdb = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM));
  
  start_time = time(NULL);
  if(strcmp(protein->parameters.algorithm, "HBUILD") != 0)
    {  

      if(strcmp(protein->parameters.algorithm, "MINIMIZE") == 0)
    	{
	  if(protein->num_ligands!=0)
	  {
	  	fprintf(stderr,"WARNING MINIMIZE w/ ligands is poorly characterized.\n");
		fprintf(stderr,"\tTherefore, ligands are fixed for now...powell can create non-orthogonal un-normalized transform matricies\n");
	  }

	  energies[0] = minimization_objective_function(delta_dihedral_array[1],&powell);
	  
	  stopflag = 0;
	  
	  j=0;
	  while(stopflag == 0)
	    {
	      ++j;
	

	      powell_minimization(delta_dihedral_array[1], num_dihedrals, protein->parameters.max_iterations, &energies[1], 
		     &minimization_objective_function,&powell, start_time);
	 

	      if(2.0*fabs(energies[0]-energies[1]) <= 2*TOL*(fabs(energies[0])+fabs(energies[1])))
		stopflag = 1;
	      else
		energies[0] = energies[1];
	      
	      now = time(NULL);
	      if(difftime(now,start_time) >= MAX_OPTIMIZATION_TIME)
		NUM_FUNCTION_CALLS = protein->parameters.max_function_calls + 2;
	      
	      if(NUM_FUNCTION_CALLS >= protein->parameters.max_function_calls)
		stopflag = 1;
	      
	      if(j >= num_dihedrals)	
		stopflag = 1;
	    }	

	  convert_delta_dihedral_array_to_CHROMOSOME(&protein->chr[2], delta_dihedral_array[1], 
												&protein->chr[1], protein->parameters.rebuild_backbone_flag, 
												protein->min_fixed_res, protein->min_float_res);
   	}
      else	/* adjust backbone torsions */
	{
	  
	  energies[0] = bkbn_torsion_overlay_rmsd(delta_dihedral_array[1],&powell);
	  
	  stopflag = 0;
	  j=0;
	  while(stopflag == 0)
	    {
	      ++j;
	      
	      powell_minimization(delta_dihedral_array[1], num_dihedrals, protein->parameters.max_iterations, &energies[1], 
		     &bkbn_torsion_overlay_rmsd,&powell, start_time);
	      if(2.0*fabs(energies[0]-energies[1]) <= 2*TOL*(fabs(energies[0])+fabs(energies[1])))
		stopflag = 1;
	      else
		energies[0] = energies[1];
	      
	      now = time(NULL);
	      if(difftime(now,start_time) >= MAX_OPTIMIZATION_TIME)
		NUM_FUNCTION_CALLS = protein->parameters.max_function_calls + 2;
	      
	      if(NUM_FUNCTION_CALLS >= protein->parameters.max_function_calls)
		stopflag = 1;
	      
	      if(j >= num_dihedrals)	
		stopflag = 1;
	    }	
	  
	  i=1;
	  protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
	  protein->chr[2].bkbngenes = protein->chr[2].first_bkbngene;
	  while(protein->chr[1].bkbngenes->nextbkbngene!=NULL)
	    {
	      protein->chr[2].bkbngenes->phi = protein->chr[1].bkbngenes->phi + delta_dihedral_array[1][i];
	      ++i;
	      if(protein->chr[2].bkbngenes->phi>180)
		protein->chr[2].bkbngenes->phi = protein->chr[2].bkbngenes->phi -360.0;
	      else
		if(protein->chr[2].bkbngenes->phi < -180)
		  protein->chr[2].bkbngenes->phi = protein->chr[2].bkbngenes->phi + 360.0;
	      
	      protein->chr[2].bkbngenes->psi = protein->chr[1].bkbngenes->psi + delta_dihedral_array[1][i];
	      ++i;
	      if(protein->chr[2].bkbngenes->psi>180)
		protein->chr[2].bkbngenes->psi = protein->chr[2].bkbngenes->psi -360.0;
	      else
		if(protein->chr[2].bkbngenes->psi < -180)
		  protein->chr[2].bkbngenes->psi = protein->chr[2].bkbngenes->psi + 360.0; 
	      
	      protein->chr[2].bkbngenes->omega = protein->chr[1].bkbngenes->omega + delta_dihedral_array[1][i];
	      ++i;
	      if(protein->chr[2].bkbngenes->omega>180)
		protein->chr[2].bkbngenes->omega = protein->chr[2].bkbngenes->omega -360.0;
	      else
		if(protein->chr[2].bkbngenes->omega < -180)
		  protein->chr[2].bkbngenes->omega = protein->chr[2].bkbngenes->omega + 360.0; 
	      
	      protein->chr[1].bkbngenes = protein->chr[1].bkbngenes->nextbkbngene;
	      protein->chr[2].bkbngenes = protein->chr[2].bkbngenes->nextbkbngene;
	    }
	  protein->chr[1].bkbngenes = protein->chr[1].first_bkbngene;
	  protein->chr[2].bkbngenes = protein->chr[2].first_bkbngene;
	  
	}

    }
    
  protein->E_working = energies[1];

  protein->final_chr = protein->chr[2];
  i = SASA_FLAG; SASA_FLAG = 1;
  final_chr_to_final_pdb_final_energy(protein); 
  SASA_FLAG = i;   

  if(protein->mini_fixed_atoms!=NULL)
    {
      free_memory(protein->mini_fixed_atoms);
      free_memory(protein->fixed_atoms);
    }
  protein->mini_fixed_atoms = mini_fixed_atoms;
  protein->fixed_atoms = fixed_atoms;
  
  if(protein->CBminus_minipdb != NULL)
    free_memory(protein->CBminus_minipdb);
  
  if(SASA_FLAG == 1 && MAX_ATOMS <= REAL_MAX_ATOMS)
    {	
      if(VECTOR_PTR!=NULL)
	{   
	  for(i=0;i<MAX_ATOMS;++i)
	    free_memory(VECTOR_PTR[i]);
	  free_memory(VECTOR_PTR);
	  VECTOR_PTR = NULL;
	}
    }	

  if(delta_dihedral_array!=NULL)
  {
  	for(i=0;i<(num_dihedrals+3);++i)
    		free_memory(delta_dihedral_array[i]);
  	free_memory(delta_dihedral_array);
  }
  if(energies!=NULL)
  	free_memory(energies);
	
}