/* EGAD: docking.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 docking one molecule to another

*/

#include "docking.h"

/* vdw, coulomb, hbond energies between arrays i_atoms and j_atoms */
/* born pair energies calc'd with born radii already in i_atoms in j_atoms */
ENERGY pair_energy_arraypair(mini_pdbATOM *i_atoms, mini_pdbATOM *j_atoms, int calc_born_flag)
{
	int i,j, vdw_contact_flag=0, contact_flag;
	VDW_ENERGY vdwE;
	double r, r2, xmax, ymax, zmax, xmin, ymin, zmin, distance_limit, distance_limit_sqrd;
	double dx, dy, dz;
	double ff_dist_cutoff;
	ENERGY energy;
	static int first_time_flag=0; 
	static ENERGY baseline_energy;
	static double *i_born=NULL, *j_born=NULL;
	extern double **VDW_CONTACT_DISTANCE, HYDROPHOBIC_ASP;

	/*
	static COULOMBIC *coulombic=NULL, *first_coulombic=NULL;
	static COULOMBIC *i_coulombic=NULL, *i_first_coulombic=NULL;
	static COULOMBIC *j_coulombic=NULL, *j_first_coulombic=NULL;
	static double **CHARGE_PRODUCT;
	double alpha, EE, f;
	*/


	/* speed things up by reducing the effective distance */
	ff_dist_cutoff = FORCEFIELD_DISTANCE_CUTOFF;
	FORCEFIELD_DISTANCE_CUTOFF = 5;
	FORCEFIELD_DISTANCE_CUTOFF_SQRD = 25;
	
	/* the first time this function is run, initialize these arrays; these will hold the born desolvation due to the
		isolated molecules stored in the "born_radius" entries of i_atoms and j_atoms 
	*/
	if(calc_born_flag == 1)
	{
		
		if(i_born==NULL)
		{
			i_born = (double *)calloc(MAX_ATOMS,sizeof(double));
			j_born = (double *)calloc(MAX_ATOMS,sizeof(double));

			i=1;
			while(i_atoms[i].seq_position!=ENDFLAG)
			{
				i_born[i] = i_atoms[i].born_radius;
				++i;
			}

			j=1;
			while(j_atoms[j].seq_position!=ENDFLAG)
			{
				j_born[j] = j_atoms[j].born_radius;
				++j;
			}

		
		/*
			if(first_coulombic == NULL)
			{	
				coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC));
				first_coulombic = coulombic;
				coulombic->next = NULL;
			}
		
			coulombic = first_coulombic;
			coulombic->r2 = ENDFLAG;


			// calculate internal distances for i_atoms and j_atoms
			i_coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC));
			i_first_coulombic = i_coulombic;
			i_coulombic->next = NULL;
			i_coulombic = i_first_coulombic;
			i_coulombic->r2 = ENDFLAG;

			j_coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC));
			j_first_coulombic = j_coulombic;
			j_coulombic->next = NULL;
			j_coulombic = j_first_coulombic;
			j_coulombic->r2 = ENDFLAG;
			
			sidechain_internal_energy(i_atoms, i_first_coulombic, i_coulombic);
			sidechain_internal_energy(j_atoms, j_first_coulombic, j_coulombic);
		*/

		
		}		
	}

	/* if the molecules are not within FORCEFIELD_DISTANCE_CUTOFF from each other, reset to 1e10 just so an
energy can be calculated and try again.....otherwise, the minimizer gets stuck. 
	This situation may come up in the early cycles where the initial random (constrained) guesses may be far off */
	
	contact_flag = 0;
	distance_limit = FORCEFIELD_DISTANCE_CUTOFF;

	while(contact_flag == 0)
	{
	distance_limit_sqrd = distance_limit*distance_limit;

	ENERGY_0(energy);
	VDW_ENERGY_0(vdwE);


	if(calc_born_flag == 1)
	{
		j=1;
		while(j_atoms[j].seq_position!=ENDFLAG)
		{
			j_atoms[j].born_radius = j_born[j];
			++j;
		}
	
		/* coulombic = first_coulombic; */

	}

	i=1;
	while(i_atoms[i].seq_position!=ENDFLAG)
	{
		if(calc_born_flag == 1)
			i_atoms[i].born_radius = i_born[i];

		if(i_atoms[i].atom_ptr->atom_ptr->atom_class!=ENDFLAG)  /* no U */
		{

		xmax=i_atoms[i].coord.x+distance_limit;
	  	xmin=i_atoms[i].coord.x-distance_limit;
	  	ymax=i_atoms[i].coord.y+distance_limit;
	  	ymin=i_atoms[i].coord.y-distance_limit;
	  	zmax=i_atoms[i].coord.z+distance_limit;
	  	zmin=i_atoms[i].coord.z-distance_limit;
	
		vdw_contact_flag=0;

		j=1;
		while(j_atoms[j].seq_position!=ENDFLAG)
		{
			if(j_atoms[j].atom_ptr->atom_ptr->atom_class!=ENDFLAG)  /* no U */
			{

		 
			if(j_atoms[j].coord.x<=xmax) 
	    		if(j_atoms[j].coord.x>=xmin)
	      		if(j_atoms[j].coord.y<=ymax)
			if(j_atoms[j].coord.y>=ymin)
		  	if(j_atoms[j].coord.z<=zmax) 
		    if(j_atoms[j].coord.z>=zmin)
			{ 
			dx = j_atoms[j].coord.x - i_atoms[i].coord.x;
			dy = j_atoms[j].coord.y - i_atoms[i].coord.y;
			dz = j_atoms[j].coord.z - i_atoms[i].coord.z;
			r2 = dx*dx + dy*dy + dz*dz; 
			if(r2>=TINY  && r2<=distance_limit_sqrd  )
			{
					contact_flag = 1;

					r=sqrt(r2);

					add_atompair_vdw_energy(&vdwE, i_atoms[i].atom_ptr, j_atoms[j].atom_ptr, r2, r);
					
			if(r<=VDW_CONTACT_DISTANCE[i_atoms[i].atom_ptr->atom_ptr->atom_class][j_atoms[j].atom_ptr->atom_ptr->atom_class])
						vdw_contact_flag=1;


					energy.E_coulomb +=  coulomb_energy_q2_per_angstrom(r,i_atoms[i].atom_ptr, j_atoms[j].atom_ptr);
					energy.E_hbond += hbond_energy(&i_atoms[i], &j_atoms[j], r, r2);
			
			 
					if(calc_born_flag == 1)
					{
						// HU does not influence the born radius 

						if(strcmp(j_atoms[j].atom_ptr->atom_ptr->atomtype, "HU")!=0)  
                                	i_atoms[i].born_radius += born_nonbond_calc(i_atoms[i].atom_ptr, j_atoms[j].atom_ptr, r2); 
                                     
 						if(strcmp(i_atoms[i].atom_ptr->atom_ptr->atomtype, "HU")!=0)  
                                	j_atoms[j].born_radius += born_nonbond_calc(j_atoms[j].atom_ptr, i_atoms[i].atom_ptr, r2); 
		
					/*
						coulombic->index_i = i;
                                coulombic->index_j = j; 
                                coulombic->r2 = r2;
                                 
						if(coulombic->next == NULL)
						{           
                                	coulombic->next = (COULOMBIC *)malloc(sizeof(COULOMBIC));
							coulombic = coulombic->next;  // mark end of linked-list 
                                	coulombic->next = NULL;
						}
						else
							coulombic = coulombic->next;

						coulombic->r2 = ENDFLAG; // mark end
					*/
				
					}
			 
	
			}
			}
			}
			++j;
		}
		}

		/*
		in lieu of SASA-based term (expensive!), penalize non-vdw-bonding hphob atoms
		if these are not in vdW contact w/ a receptor atom, they are probably not buried, 
		and they don't bury a receptor atom (hence the multiply by 2...assume that it would bury an equal area on the receptor)
		probably an overestimate, but a reasonable estimate nonetheless
		the sasa here is the sasa of the free state
		*/

		if(vdw_contact_flag==0)
			if(i_atoms[i].atom_ptr->atom_ptr->asp > 0)
				energy.E_sasa += i_atoms[i].sasa;

		if(calc_born_flag == 1)
		{
			if(i_atoms[i].born_radius >= 0)
                       i_atoms[i].born_radius = -1.0;
                    i_atoms[i].born_radius = -HALF_COULOMB_CONST/i_atoms[i].born_radius;

                	energy.E_born += BORN_CONST(i_atoms[i].born_radius)*i_atoms[i].atom_ptr->charge2/i_atoms[i].born_radius; 

		}
		
		++i;
	}

	if(contact_flag==0)
		distance_limit = 1e10;
		

	}


	if(calc_born_flag == 1)
	{
		j=1;
		while(j_atoms[j].seq_position!=ENDFLAG)
		{
			if(j_atoms[j].born_radius >= 0)
                       j_atoms[j].born_radius = -1.0;
                    j_atoms[j].born_radius = -HALF_COULOMB_CONST/j_atoms[j].born_radius;

                	energy.E_born += BORN_CONST(j_atoms[j].born_radius)*j_atoms[j].atom_ptr->charge2/j_atoms[j].born_radius; 

			++j;
		}

	/*
		coulombic = first_coulombic;
		while(coulombic->r2!=ENDFLAG)
		{
			alpha = i_atoms[coulombic->index_i].born_radius*j_atoms[coulombic->index_j].born_radius;
			f = f_of_r2_alpha(coulombic->r2, alpha);
	EE = CHARGE_PRODUCT[i_atoms[coulombic->index_i].atom_ptr->charge_class][j_atoms[coulombic->index_j].atom_ptr->charge_class]/f; 

			if(EE > POL_ATTRACTION_CAP)        
                		EE = POL_ATTRACTION_CAP;
							
			energy.E_pol += BORN_CONST(f)*EE;

			coulombic = coulombic->next;
		}

		i_coulombic = i_first_coulombic;
		while(i_coulombic->r2!=ENDFLAG)
		{
			alpha = i_atoms[i_coulombic->index_i].born_radius*i_atoms[i_coulombic->index_j].born_radius;
			f = f_of_r2_alpha(i_coulombic->r2, alpha);
	EE = CHARGE_PRODUCT[i_atoms[i_coulombic->index_i].atom_ptr->charge_class][i_atoms[i_coulombic->index_j].atom_ptr->charge_class]/f; 

			if(EE > POL_ATTRACTION_CAP)        
                		EE = POL_ATTRACTION_CAP;
							
			energy.E_pol += BORN_CONST(f)*EE;

			i_coulombic = i_coulombic->next;
		}

		j_coulombic = j_first_coulombic;
		while(j_coulombic->r2!=ENDFLAG)
		{
			alpha = j_atoms[j_coulombic->index_i].born_radius*j_atoms[j_coulombic->index_j].born_radius;
			f = f_of_r2_alpha(j_coulombic->r2, alpha);
	EE = CHARGE_PRODUCT[j_atoms[j_coulombic->index_i].atom_ptr->charge_class][j_atoms[j_coulombic->index_j].atom_ptr->charge_class]/f; 

			if(EE > POL_ATTRACTION_CAP)        
                		EE = POL_ATTRACTION_CAP;
							
			energy.E_pol += BORN_CONST(f)*EE;

			j_coulombic = j_coulombic->next;
		}
	*/
			
	}

	energy.E_vdw = sum_VDW_ENERGY(vdwE);
	energy.clashes = vdwE.clashes;
	energy.E_coulomb = COULOMB_CONST_OVER_INTERNAL_DIELECTRIC*energy.E_coulomb;
	energy.E_pol = COULOMB_CONST*energy.E_pol;   
        energy.E_born = HALF_COULOMB_CONST*energy.E_born;
	energy.E_sasa = 2*HYDROPHOBIC_ASP*energy.E_sasa;   

	if(first_time_flag==0)
	{
		baseline_energy = energy;
		first_time_flag = 1;
	}

	energy = ENERGY_SUBTRACT(energy, baseline_energy);

	FORCEFIELD_DISTANCE_CUTOFF = ff_dist_cutoff;
	FORCEFIELD_DISTANCE_CUTOFF_SQRD = FORCEFIELD_DISTANCE_CUTOFF*FORCEFIELD_DISTANCE_CUTOFF;


	return(energy);
}

double rigid_dock_objective_function(double *translate_rotate_array, POWELL *powell)
{
	int i;
	ENERGY energy;
	static double *transform_matrix=NULL;
	static mini_pdbATOM *transformed_pdb=NULL;
	static int max_atoms=0;
	extern int NUM_FUNCTION_CALLS, GB_FLAG;
	

	/* check bounds .... penalize severely if outside */
	for(i=1;i<=7;++i)
	{
		if(translate_rotate_array[i] > powell->rigid_docking_data->upper_bound[i])
			return(1e10);
		if(translate_rotate_array[i] < powell->rigid_docking_data->lower_bound[i])
			return(1e10);
	}

	
	
	++NUM_FUNCTION_CALLS;

	if(transform_matrix == NULL)
		transform_matrix = (double *)calloc(14,sizeof(double));
		
	if(max_atoms<MAX_ATOMS)
	{
		max_atoms = MAX_ATOMS;
		if(transformed_pdb!=NULL)
			free_memory(transformed_pdb);

		transformed_pdb = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));

		i=1;
		while(powell->rigid_docking_data->starting_pdb[i].seq_position!=ENDFLAG)
		{
			transformed_pdb[i] = powell->rigid_docking_data->starting_pdb[i];
			++i;
		}
		transformed_pdb[i] = powell->rigid_docking_data->starting_pdb[i];

	}


	
	translate_rotate_array_to_transform_matrix(translate_rotate_array, transform_matrix);

	i=1;
	while(powell->rigid_docking_data->starting_pdb[i].seq_position!=ENDFLAG)
		{
			transformed_pdb[i].coord = tranform_coordinates(powell->rigid_docking_data->starting_pdb[i].coord,transform_matrix);
			++i;
		}
	
	/* use GB only if all the atoms are real; if var_pos not NULL, then pseudoatoms are included */
	if(powell->protein->var_pos==NULL)
		energy = pair_energy_arraypair(transformed_pdb, powell->rigid_docking_data->docking_target_pdb, GB_FLAG);
	else
		energy = pair_energy_arraypair(transformed_pdb, powell->rigid_docking_data->docking_target_pdb, 0);

	

	energy.E_total = ENERGY_TOTAL(energy);

	if(energy.E_total < powell->rigid_docking_data->energy)
	{
		powell->rigid_docking_data->energy = energy.E_total;
		for(i=1;i<=7;++i)
				powell->rigid_docking_data->translate_rotate_array[i] = translate_rotate_array[i];


		// printf("%d\t%lf\n",NUM_FUNCTION_CALLS,powell->rigid_docking_data->energy);
	}

	return(energy.E_total);

}





void rigid_dock(PROTEIN *protein)
{
	div_t chain_id_number;
    	char chain_id;
	time_t start_time;
	mini_pdbATOM *free_A, *free_B, *initial_complex_pdb, *target_pdb;
    	int i,a,b,j,n,found_flag,k,A_flag, uul_index, psl_index, atomname_index,found_it_flag;
	int moving_start_seqpos, moving_end_seqpos;
	double dihed_CB, r_max, r;
	char *chain_id_list;
	SUPER_CHAIN_ID_LIST *super_chain_id_list;
	POWELL *powell;
	CARTESIAN centroid;
	double min_theta, min_dx, min_dy, min_dz, max_theta, max_dx, max_dy, max_dz;
	VARIABLE_POSITION *dummy_varpos;
	extern int NUM_FUNCTION_CALLS;
	extern ATOMRESPARAM *U_ATOMRESPARAM;

	
	uul_index = search_residuetype_RESPARAM("UUL", protein->resparam);
	psl_index = search_atomname_ATOMRESPARAM("PSL", protein->resparam[uul_index].sorted_atom);


	initial_complex_pdb = (mini_pdbATOM *)calloc(MAX_ATOMS,sizeof(mini_pdbATOM));
	i=1;
	while(protein->mini_Template[i].seq_position != ENDFLAG)
	{
		initial_complex_pdb[i] = protein->mini_Template[i];
		++i;
	}
	initial_complex_pdb[i] = protein->mini_Template[i];

	

	//  if varpos is defined replace moving sidechains (not ligands) for varpos with UUL 
	if(protein->var_pos!=NULL)
	{
		i=1; k=1; j=1; n=1;
		while(protein->var_pos[i].seq_position!=ENDFLAG)
		{
			if(protein->var_pos[i].fixed_flag == 0 && protein->var_pos[i].choice[1].resparam_ptr->ligand_flag == 0)
			{
				while(protein->var_pos[i].seq_position != protein->mini_Template[n].seq_position)
				{
					initial_complex_pdb[k] = protein->mini_Template[n];
					++k; ++n;
				}	
				
				while(protein->var_pos[i].seq_position == protein->mini_Template[n].seq_position)
				{
					atomname_index = search_atomname_ATOMRESPARAM(protein->mini_Template[n].atom_ptr->atomname, 
																	protein->resparam[uul_index].sorted_atom);

					if(atomname_index != 0)
					{
						initial_complex_pdb[k] = protein->mini_Template[n];
						initial_complex_pdb[k].atom_ptr = &(protein->resparam[uul_index].sorted_atom[atomname_index]);
						++k; 
					}					

					++n;
				}

				/* append PSL */
				initial_complex_pdb[k] = initial_complex_pdb[k-1]; 
				initial_complex_pdb[k].atom_ptr = &(protein->resparam[uul_index].sorted_atom[psl_index]);

				while(protein->bkbn[j].seq_position!=protein->var_pos[i].seq_position)
					++j;
	
				dihed_CB = protein->bkbn[j].phi-120;
				initial_complex_pdb[k].coord = chi2xyz(&protein->bkbn[j].CA, &protein->bkbn[j].N, &protein->bkbn[j-1].C, 
                                                  &initial_complex_pdb[k].atom_ptr->bondlength, 
                                                  &initial_complex_pdb[k].atom_ptr->bondangle, 
                                                  &dihed_CB);
				++k;
				 

			}
			++i;
		}
		while(protein->mini_Template[n].seq_position!=ENDFLAG)
		{
			initial_complex_pdb[k] = protein->mini_Template[n];
			++k; ++n;
		}
		initial_complex_pdb[k].atom_ptr = NULL;
		initial_complex_pdb[k].seq_position=ENDFLAG;
	}
	
	powell = (POWELL *)malloc(sizeof(POWELL));
	powell->protein = protein;
	
	powell->rigid_docking_data = (RIGID_DOCKING_DATATYPE *)malloc(sizeof(RIGID_DOCKING_DATATYPE));
	
	free_A = (mini_pdbATOM *)calloc(MAX_ATOMS,sizeof(mini_pdbATOM));
	free_B  = (mini_pdbATOM *)calloc(MAX_ATOMS,sizeof(mini_pdbATOM));

	// find the smaller subunit....this will be moved; the larger will remain stationary 

	chain_id_list = protein->chain_id_list;
	super_chain_id_list = protein->super_chain_id_list;


  	i=1; k=0; a=1;b=1;
        while(initial_complex_pdb[i].seq_position != ENDFLAG)
        {
        	chain_id_number = div(initial_complex_pdb[i].seq_position, 1000);
              	chain_id = chain_id_list[chain_id_number.quot+1];

                found_flag=0; k=1;
                while(super_chain_id_list[1].chain_id[k]!='\0' && found_flag==0)
                {
			if(chain_id == super_chain_id_list[1].chain_id[k])
				found_flag=1;
			++k;
                }

		if(found_flag==1)       // count the two superchains 
			++a;
		else
			++b;

		++i;
         }

         if(a<b)         // if the defined superchain is the smaller one, copy into free_A
               A_flag=1;
         else
               A_flag=0;       // otherwise, put the other into free_A


	split_complex(initial_complex_pdb, free_A, free_B, chain_id_list, super_chain_id_list);

	if(A_flag==1)
	{
		powell->rigid_docking_data->starting_pdb = free_A;
		powell->rigid_docking_data->docking_target_pdb = free_B;
	}
        else
	{
		powell->rigid_docking_data->starting_pdb = free_B;
		powell->rigid_docking_data->docking_target_pdb = free_A;
	}
	
	/* translate centroid of starting pdb to the origin...makes interpreting rotations much easier */
	centroid = get_centroid(powell->rigid_docking_data->starting_pdb);
	centroid.x = -centroid.x;
	centroid.y = -centroid.y;
	centroid.z = -centroid.z;
	i=1;
	while(powell->rigid_docking_data->starting_pdb[i].seq_position!=ENDFLAG)
	{
		powell->rigid_docking_data->starting_pdb[i].coord = 
				add_CARTESIAN(powell->rigid_docking_data->starting_pdb[i].coord, centroid);
		++i;
	}
	centroid.x = -centroid.x;
	centroid.y = -centroid.y;
	centroid.z = -centroid.z;


	/* calc SASA and born desolvation factors */
	energy_calc(powell->rigid_docking_data->starting_pdb, 3, 1);
	energy_calc(powell->rigid_docking_data->docking_target_pdb, 3, 1);

	moving_start_seqpos =  powell->rigid_docking_data->starting_pdb[1].seq_position;
	i=1;
	while(powell->rigid_docking_data->starting_pdb[i].seq_position!=ENDFLAG)
		++i;
	moving_end_seqpos = powell->rigid_docking_data->starting_pdb[i-1].seq_position;


	// minimize the initial docked structure

		min_theta = -18; max_theta = 18;

		min_dx = -5 + centroid.x; max_dx = 5 + centroid.x; 
		min_dy = -5 + centroid.y; max_dy = 5 + centroid.y; 
		min_dz = -5 + centroid.z; max_dz = 5 + centroid.z; 
		
	
	if(protein->parameters.dock_local_flag == 0)	// de novo docking
	{
		min_theta = -180; max_theta = 180;

		target_pdb = (mini_pdbATOM *)calloc(MAX_ATOMS,sizeof(mini_pdbATOM));

		if(protein->var_pos!=NULL)	// keep near variable positions
		{
			i=1; k=1; n=1;
			while(protein->var_pos[i].seq_position!=ENDFLAG)
			{
				if(protein->var_pos[i].fixed_flag == 0 && protein->var_pos[i].choice[1].resparam_ptr->ligand_flag == 0)
				{
					while(powell->rigid_docking_data->docking_target_pdb[n].seq_position != protein->var_pos[i].seq_position &&
							powell->rigid_docking_data->docking_target_pdb[n].seq_position!=ENDFLAG)
						++n;

					if(powell->rigid_docking_data->docking_target_pdb[n].seq_position!=ENDFLAG)
					{
					while(powell->rigid_docking_data->docking_target_pdb[n].seq_position == protein->var_pos[i].seq_position)
					{
						target_pdb[k] = powell->rigid_docking_data->docking_target_pdb[n];
						++n; ++k;
					}
					}
				}
				++i;
			}
			target_pdb[k].seq_position=ENDFLAG; target_pdb[k].atom_ptr=NULL;

		}
		else	// anything goes
		{
			n=1;
			while(powell->rigid_docking_data->docking_target_pdb[n].seq_position != ENDFLAG)
			{
				target_pdb[n] = powell->rigid_docking_data->docking_target_pdb[n];
				++n;
			}
			target_pdb[n] = powell->rigid_docking_data->docking_target_pdb[n];

		}


		// define min and max translations
		// define cube that surrounds the receptor plus 1/2 ligand-length
		// define max displacements 
		
		// find length of moving molec
		r_max = 0;
		i=1;
		while(powell->rigid_docking_data->starting_pdb[i].seq_position!=ENDFLAG)
		{
			if(powell->rigid_docking_data->starting_pdb[i].atom_ptr->atom_ptr->atom_class!=ENDFLAG)
			{
			j=i+1;
			while(powell->rigid_docking_data->starting_pdb[j].seq_position!=ENDFLAG)
			{
				if(powell->rigid_docking_data->starting_pdb[j].atom_ptr->atom_ptr->atom_class!=ENDFLAG)
				{
					r = 
	Distance_sqrd(powell->rigid_docking_data->starting_pdb[i].coord, powell->rigid_docking_data->starting_pdb[j].coord);
			
					if(r > r_max)
						r_max = r;
				}
				++j;
			}
			}
			++i;
		}

		r_max = sqrt(r_max)/2.0;

		// find min_xyz and max_xyz on target

		min_dx = 1e10; max_dx = -1e10; 
		min_dy = 1e10; max_dy = -1e10; 
		min_dz = 1e10; max_dz = -1e10; 

		i=1;
		while(target_pdb[i].seq_position!=ENDFLAG)
		{
			if(min_dx > target_pdb[i].coord.x)
				min_dx = target_pdb[i].coord.x;
			if(min_dy > target_pdb[i].coord.y)
				min_dy = target_pdb[i].coord.y;
			if(min_dz > target_pdb[i].coord.z)
				min_dz = target_pdb[i].coord.z;
						

			if(max_dx < target_pdb[i].coord.x)
				max_dx = target_pdb[i].coord.x;
			if(max_dy < target_pdb[i].coord.y)
				max_dy = target_pdb[i].coord.y;
			if(max_dz < target_pdb[i].coord.z)
				max_dz = target_pdb[i].coord.z;
			

			++i;
		}
		
	 
		min_dx =  min_dx - r_max;
		min_dy =  min_dy - r_max;
		min_dz =  min_dz - r_max;

		max_dx = max_dx + r_max;
		max_dy = max_dy + r_max;
		max_dz = max_dz + r_max;
	 
		 /*
		printf("%lf\t%lf\t%lf\n",min_dx,min_dy,min_dz);
		printf("%lf\t%lf\t%lf\n\n",max_dx,max_dy,max_dz);
		 */

		if(min_dx < max_dx)
		{
			min_dx = min_dx;
			max_dx = max_dx;
		}
		else
		{
			r = min_dx;
			min_dx = max_dx;
			max_dx = r;
		}

		if(min_dy  < max_dy )
		{
			min_dy = min_dy;
			max_dy = max_dy;
		}
		else
		{
			r = min_dy;
			min_dy = max_dy;
			max_dy = r;
		}

		if(min_dz  < max_dz )
		{
			min_dz = min_dz;
			max_dz = max_dz;
		}
		else
		{
			r = min_dz;
			min_dz = max_dz;
			max_dz = r;
		}


		 /*
		printf("%lf\t%lf\t%lf\n",min_dx,min_dy,min_dz);
		printf("%lf\t%lf\t%lf\n",max_dx,max_dy,max_dz);
		 */

		free_memory(target_pdb);

	}
	if(protein->parameters.dock_local_flag == 2) // minimize the initial docked structure, but allow large rotations
	{
		min_theta = -180; max_theta = 180;
		min_dx = -10 + centroid.x; max_dx = 10 + centroid.x; 
		min_dy = -10 + centroid.y; max_dy = 10 + centroid.y; 
		min_dz = -10 + centroid.z; max_dz = 10 + centroid.z; 
	}

	
	powell->rigid_docking_data->translate_rotate_array = (double *)calloc(9,sizeof(double));
	
	powell->rigid_docking_data->lower_bound = (double *)calloc(9,sizeof(double));
	powell->rigid_docking_data->upper_bound = (double *)calloc(9,sizeof(double));

	powell->rigid_docking_data->translate_rotate_array[1] = centroid.x;	
		powell->rigid_docking_data->lower_bound[1] = min_dx; powell->rigid_docking_data->upper_bound[1] = max_dx; // dx
	powell->rigid_docking_data->translate_rotate_array[2] = centroid.y;	
		powell->rigid_docking_data->lower_bound[2] = min_dy; powell->rigid_docking_data->upper_bound[2] = max_dy; // dy
	powell->rigid_docking_data->translate_rotate_array[3] = centroid.z;	
		powell->rigid_docking_data->lower_bound[3] = min_dz; powell->rigid_docking_data->upper_bound[3] = max_dz; // dz
	
	powell->rigid_docking_data->translate_rotate_array[4] = 0.1;	
		powell->rigid_docking_data->lower_bound[4] = -1; powell->rigid_docking_data->upper_bound[4] = 1; // i
	powell->rigid_docking_data->translate_rotate_array[5] = 0.1;	
		powell->rigid_docking_data->lower_bound[5] = -1; powell->rigid_docking_data->upper_bound[5] = 1; // j
	powell->rigid_docking_data->translate_rotate_array[6] = 0.1;	
		powell->rigid_docking_data->lower_bound[6] = -1; powell->rigid_docking_data->upper_bound[6] = 1; // k

	// rotation angle around <i,j,k>
	powell->rigid_docking_data->translate_rotate_array[7] = 0;	
		powell->rigid_docking_data->lower_bound[7] = min_theta; powell->rigid_docking_data->upper_bound[7] = max_theta; 

	
	powell->rigid_docking_data->energy = rigid_dock_objective_function(powell->rigid_docking_data->translate_rotate_array, powell);
	protein->E_working = powell->rigid_docking_data->energy;

	NUM_FUNCTION_CALLS = 0;
	powell->protein->parameters.max_function_calls = (int)2e6;
	start_time = time(NULL);
	
		general_purpose_GA(powell->rigid_docking_data->translate_rotate_array, powell->rigid_docking_data->lower_bound, powell->rigid_docking_data->upper_bound, 7, protein->parameters.ga_convergence, &protein->E_working, &rigid_dock_objective_function, powell, start_time);		

	/*
		powell_minimization(powell->rigid_docking_data->translate_rotate_array, 7, protein->parameters.max_iterations, &protein->E_working, &rigid_dock_objective_function, powell, start_time);
	 */

		

	protein->E_working = powell->rigid_docking_data->energy;

	if(protein->transform_matrix==NULL)
		protein->transform_matrix = (double *)calloc(14,sizeof(double));

	translate_rotate_array_to_transform_matrix(powell->rigid_docking_data->translate_rotate_array, protein->transform_matrix);
	
	
	

	// transform the moved chain
	
	j=1;
	while(powell->rigid_docking_data->starting_pdb[j].seq_position!=ENDFLAG)
	{

		i=1;	
		while(powell->rigid_docking_data->starting_pdb[j].seq_position != protein->Template[i].seq_position)
			++i;

			/* loop through the residue and generate coords for each atom or set to U */
			n=i;
			while(protein->Template[n].seq_position == protein->Template[i].seq_position)
			{
				k=j; found_it_flag=0;
				while(powell->rigid_docking_data->starting_pdb[j].seq_position == 
						powell->rigid_docking_data->starting_pdb[k].seq_position)
				{

	if(strcmp(powell->rigid_docking_data->starting_pdb[k].atom_ptr->atomname, protein->Template[n].atom_ptr->atomname)==0)
					{
						protein->Template[n].coord = 
			tranform_coordinates(powell->rigid_docking_data->starting_pdb[k].coord,protein->transform_matrix);
						found_it_flag=1;
					}

					++k;
				}	

				k=j;				

				if(found_it_flag == 0)
				{
					
					strcpy(protein->Template[n].atomname,"U");
					protein->Template[n].atom_ptr = U_ATOMRESPARAM;

				}

				++n;
			}

		while(powell->rigid_docking_data->starting_pdb[j].seq_position == 
						powell->rigid_docking_data->starting_pdb[k].seq_position)
			++j;
	}

	pdbATOM_to_CHROMOSOME(protein->Template, &protein->final_chr, protein->resparam);

	
	pdbATOM_to_mini_pdbATOM(protein->Template,protein->mini_Template);

	dummy_varpos = (VARIABLE_POSITION *)calloc(MAX_RESIDUES, sizeof(VARIABLE_POSITION));
	protein->fixed_atoms = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM));
	
            i=1; j=0;
            while(protein->Template[i].seq_position!=ENDFLAG)
            {
                ++j;
                dummy_varpos[j].seq_position = protein->Template[i].seq_position;
                while(protein->Template[i].seq_position == dummy_varpos[j].seq_position)
                    ++i;

            }
            dummy_varpos[j+1].seq_position = ENDFLAG;

            make_fixed_pdbATOM(protein->Template, dummy_varpos, protein->fixed_atoms);

	protein->mini_fixed_atoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
	pdbATOM_to_mini_pdbATOM(protein->fixed_atoms,protein->mini_fixed_atoms);


	final_chr_to_final_pdb_final_energy(protein);


	free_memory(powell->rigid_docking_data->lower_bound); free_memory(powell->rigid_docking_data->upper_bound); 
	free_memory(powell->rigid_docking_data->translate_rotate_array);
	free_memory(powell->rigid_docking_data);
	free_memory(powell);


	free_memory(free_B); 
	free_memory(free_A);  
	free_memory(initial_complex_pdb); 
 

}


double rmsd_align_objective_function(double *translate_rotate_array, POWELL *powell)
{
	extern int NUM_FUNCTION_CALLS;
	int i;
	double rmsd;
	static double *transform_matrix=NULL;
	CARTESIAN transformed_coord;
	
	++NUM_FUNCTION_CALLS;
	
	if(transform_matrix==NULL)
		transform_matrix = (double *)calloc(14,sizeof(double));	
	
	translate_rotate_array_to_transform_matrix(translate_rotate_array, transform_matrix);
	
	rmsd = 0;
	for(i=1;i<=powell->rmsd_overlay_data->num_points;++i)
	{
		transformed_coord = tranform_coordinates(powell->rmsd_overlay_data->initial_coord[i], transform_matrix);
		rmsd += Distance_sqrd(transformed_coord,powell->rmsd_overlay_data->target_coord[i]);
	}
	
	rmsd = sqrt(rmsd/powell->rmsd_overlay_data->num_points);

	if(rmsd<powell->rmsd_overlay_data->rmsd)
	{
		powell->rmsd_overlay_data->rmsd = rmsd;
		
		for(i=1;i<=7;++i)
			powell->rmsd_overlay_data->translate_rotate_array[i] = translate_rotate_array[i];
	}
	
	return(rmsd);
}




void parse_pdbline(char *line, char *restype, char *atomname, char *seqpos_text, CARTESIAN *coord)
{
	int a,b, seq_pos;
	char chain_id;
	static char *dummystring=NULL;	

	if(dummystring == NULL)
		dummystring = (char *)calloc(30,sizeof(char));


	for(a=0;a<20;++a)
                dummystring[a] = '\0';
             	a=0; b = 12;
             	while(b<=15)
             	{   
                	dummystring[a] =  line[b];
                	++a; ++b; 
             	}
             	sscanf(dummystring, "%s", atomname);
		
		for(a=0;a<20;++a)
                	dummystring[a] = '\0';
             	a=0; b = 17;
             	while(b<=19)
             	{   
                	dummystring[a] =  line[b];
                	++a; ++b; 
             	}
             	sscanf(dummystring, "%s", restype);

		for(a=0;a<20;++a)
                	dummystring[a] = '\0';
             	a=0; b = 22;
             	while(b<=25)
             	{   
                	dummystring[a] =  line[b];
                	++a; ++b; 
             	}
             	sscanf(dummystring, "%d", &seq_pos);
             	
		for(a=0;a<20;++a)
                	dummystring[a] = '\0';
             	a=0; b = 21;
             	while(b<=21)
             	{   
                	dummystring[a] =  line[b];
                	++a; ++b; 
             	}
             	sscanf(dummystring, "%c", &chain_id);
		
		if(chain_id!=' ')
			sprintf(seqpos_text,"%d%c",seq_pos,chain_id);
		else
			sprintf(seqpos_text,"%d",seq_pos);

		for(a=0;a<20;++a)
                	dummystring[a] = '\0';
             	a=0; b = 30;
             	while(b<=37)
             	{   
                	dummystring[a] =  line[b];
                	++a; ++b; 
             	}
                sscanf(dummystring, "%lf", &coord->x);
             
             	for(a=0;a<20;++a)
                	dummystring[a] = '\0';
             	a=0; b = 38;
             	while(b<=45)
             	{   
                	dummystring[a] =  line[b];
                	++a; ++b; 
             	}
                sscanf(dummystring, "%lf", &coord->y);
		
		 for(a=0;a<20;++a)
                	dummystring[a] = '\0';
             	a=0; b = 46;
             	while(b<=53)
             	{   
               		dummystring[a] =  line[b];
                	++a; ++b; 
             	}
                sscanf(dummystring, "%lf", &coord->z);

}

void align_structures(PROTEIN *protein)
{
	int i,j, iteration_ctr;	
	FILE *file_ptr;
	long fp;
	double previous_rmsd;
	char *line, *restype, *atomname, *seqpos_text;
	POWELL *powell;
	VARIABLE_POSITION *dummy_varpos;
	time_t start_time;
	extern int NUM_FUNCTION_CALLS, OUTPUT_ENERGY_PER_ATOM_FLAG;

	line = (char *)calloc(MAXLINE,sizeof(char));
	restype = (char *)calloc(30,sizeof(char)); 
	atomname = (char *)calloc(30,sizeof(char));
	seqpos_text = (char *)calloc(30,sizeof(char));
	
	OUTPUT_ENERGY_PER_ATOM_FLAG = 0;
	
	if(protein->first_align_def == NULL)
		failure_report("ALIGNMENT_DEFINITION not defined for JOBTYPE ALIGN_STRUCTURES","exit");
	
	powell = (POWELL *)malloc(sizeof(POWELL));

	powell->rmsd_overlay_data = (RMSD_OVERLAY_DATATYPE *)malloc(sizeof(RMSD_OVERLAY_DATATYPE));
	
	powell->rmsd_overlay_data->num_points=0;
	protein->align_def = protein->first_align_def;
	while(protein->align_def->next!=NULL)
	{
		++powell->rmsd_overlay_data->num_points;
		protein->align_def = protein->align_def->next;
	}
	protein->num_align_points = powell->rmsd_overlay_data->num_points;

	if(protein->num_align_points < 4)
		failure_report("ERROR JOBTYPE ALIGN_STRUCTURES requires at least 4 defined overlay points","exit");

	powell->protein = protein;

	powell->rmsd_overlay_data->initial_coord = 
		(CARTESIAN *)calloc(powell->rmsd_overlay_data->num_points +2 ,sizeof(CARTESIAN));

	powell->rmsd_overlay_data->target_coord = 
		(CARTESIAN *)calloc(powell->rmsd_overlay_data->num_points +2 ,sizeof(CARTESIAN));

	

	/* open alignment target pdb file and advance to the ATOM section */
	file_ptr = fopen_file(protein->competitor_Template[1],"r");
	fp = ftell(file_ptr);
	fgets(line,MAXLINE,file_ptr);
	while(strncmp(line,"ATOM",4)!=0 && strncmp(line,"LIG",3)!=0)
	{
		fp = ftell(file_ptr);
		fgets(line,MAXLINE,file_ptr);
	}
	
	
	protein->align_def = protein->first_align_def;
	i=1;
	while(protein->align_def->next!=NULL)
	{
		fseek(file_ptr, fp, 0);	/* reset file ptr to the start of the ATOM section */
	
		fgets(line,MAXLINE,file_ptr);
		parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->target_coord[i]);


		// advance to the residue of interest

		if(strstr(protein->align_def->seqpos_text_target,"l")==0)  // protein atom target
		{
		while(strcmp(seqpos_text,protein->align_def->seqpos_text_target)!=0)
		{
			if(fgets(line,MAXLINE,file_ptr)==NULL)
			{
				sprintf(line,"ERROR Cannot find position %s in %s for alignment\n",protein->align_def->seqpos_text_target,protein->competitor_Template[1]);
				failure_report(line,"exit");
			}
			if(strncmp(line,"END",3)==0)
			{
				sprintf(line,"ERROR Cannot find position %s in %s for alignment\n",protein->align_def->seqpos_text_target,protein->competitor_Template[1]);
				failure_report(line,"exit");
			}		

			parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->target_coord[i]);
		}
		}
		else // ligand atom target....advance to the ligand
		{
		while(strcmp(restype,protein->align_def->residuetype_target)!=0)
		{
			if(fgets(line,MAXLINE,file_ptr)==NULL)
			{
				sprintf(line,"ERROR Cannot find ligand %s in %s for alignment\n",protein->align_def->residuetype_target,protein->competitor_Template[1]);
				failure_report(line,"exit");
			}
			if(strncmp(line,"END",3)==0)
			{
				sprintf(line,"ERROR Cannot find ligand %s in %s for alignment\n",protein->align_def->residuetype_target,protein->competitor_Template[1]);
				failure_report(line,"exit");
			}		

			parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->target_coord[i]);
		}
		}
		while(strcmp(atomname,protein->align_def->atomname_target)!=0)  // scan through this residue to find the atom of interest
		{
			if(fgets(line,MAXLINE,file_ptr)==NULL)
			{
				sprintf(line,"ERROR Cannot find atom %s for position %s in %s for alignment\n",protein->align_def->atomname_target, protein->align_def->seqpos_text_target,protein->competitor_Template[1]);
				failure_report(line,"exit");
			}
			if(strncmp(line,"END",3)==0)
			{
				sprintf(line,"ERROR Cannot find atom %s for position %s in %s for alignment\n",protein->align_def->atomname_target, protein->align_def->seqpos_text_target,protein->competitor_Template[1]);
				failure_report(line,"exit");
			}
			
			parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->target_coord[i]);
		}
		
		++i;
		protein->align_def = protein->align_def->next;
	}
	protein->align_def = protein->first_align_def;

	fclose(file_ptr);

	
	/* open  pdb file and advance to the ATOM section */
	file_ptr = fopen_file(protein->Template_filename,"r");
	
	fp = ftell(file_ptr);
	fgets(line,MAXLINE,file_ptr);
	while(strncmp(line,"ATOM",4)!=0 && strncmp(line,"LIG",3)!=0)
	{
		fgets(line,MAXLINE,file_ptr);
		fp = ftell(file_ptr);
	}	

	protein->align_def = protein->first_align_def;
	i=1;
	while(protein->align_def->next!=NULL)
	{
		fseek(file_ptr, fp, 0);	/* reset file ptr to the start of the ATOM section */
	
		fgets(line,MAXLINE,file_ptr);
		parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->initial_coord[i]);

		if(strstr(protein->align_def->seqpos_text_initial,"l")==0)
		{
		while(strcmp(seqpos_text,protein->align_def->seqpos_text_initial)!=0)
		{
			if(fgets(line,MAXLINE,file_ptr)==NULL)
			{
				sprintf(line,"ERROR Cannot find position %s in %s for alignment\n",protein->align_def->seqpos_text_initial,protein->Template_filename);
				failure_report(line,"exit");
			}
			if(strncmp(line,"END",3)==0)
			{
				sprintf(line,"ERROR Cannot find position %s in %s for alignment\n",protein->align_def->seqpos_text_initial,protein->Template_filename);
				failure_report(line,"exit");
			}		

			parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->initial_coord[i]);
		}
		}
		else  // for ligands, advance to the ligand
		{
		while(strcmp(restype,protein->align_def->residuetype_initial)!=0)
		{
			if(fgets(line,MAXLINE,file_ptr)==NULL)
			{
				sprintf(line,"ERROR Cannot find ligand %s in %s for alignment\n",protein->align_def->residuetype_initial,protein->Template_filename);
				failure_report(line,"exit");
			}
			if(strncmp(line,"END",3)==0)
			{
				sprintf(line,"ERROR Cannot find ligand %s in %s for alignment\n",protein->align_def->residuetype_initial,protein->Template_filename);
				failure_report(line,"exit");
			}		

			parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->initial_coord[i]);
		}
		}
		while(strcmp(atomname,protein->align_def->atomname_initial)!=0)
		{
			if(fgets(line,MAXLINE,file_ptr)==NULL)
			{
				sprintf(line,"ERROR Cannot find atom %s for position %s in %s for alignment\n",protein->align_def->atomname_initial, protein->align_def->seqpos_text_initial,protein->Template_filename);
				failure_report(line,"exit");
			}
			if(strncmp(line,"END",3)==0)
			{
				sprintf(line,"ERROR Cannot find atom %s for position %s in %s for alignment\n",protein->align_def->atomname_initial, protein->align_def->seqpos_text_initial,protein->Template_filename);
				failure_report(line,"exit");
			}

			parse_pdbline(line, restype, atomname, seqpos_text, &powell->rmsd_overlay_data->initial_coord[i]);
		}


		++i;
		protein->align_def = protein->align_def->next;
	}
	protein->align_def = protein->first_align_def;
	fclose(file_ptr);

	free_memory(restype); free_memory(atomname); free_memory(seqpos_text);
	free_memory(line);

	powell->rmsd_overlay_data->translate_rotate_array = (double *)calloc(9,sizeof(double));

	powell->rmsd_overlay_data->translate_rotate_array[1] = 0;	// dx
	powell->rmsd_overlay_data->translate_rotate_array[2] = 0;	// dy
	powell->rmsd_overlay_data->translate_rotate_array[3] = 0;	// dz
	
	powell->rmsd_overlay_data->translate_rotate_array[4] = 1;	// i
	powell->rmsd_overlay_data->translate_rotate_array[5] = 1;	// j
	powell->rmsd_overlay_data->translate_rotate_array[6] = 1;	// k
	powell->rmsd_overlay_data->translate_rotate_array[7] = 0;	// rotation angle around <i,j,k>

	powell->rmsd_overlay_data->rmsd = rmsd_align_objective_function(powell->rmsd_overlay_data->translate_rotate_array, powell);
	protein->align_rmsd = powell->rmsd_overlay_data->rmsd;

	NUM_FUNCTION_CALLS = 0;
	protein->parameters.max_iterations = 1000;
	powell->protein->parameters.max_function_calls = 1000*protein->parameters.max_iterations;
	start_time = time(NULL);
	previous_rmsd = 1e10;

	iteration_ctr=0;
	while(NUM_FUNCTION_CALLS<=powell->protein->parameters.max_function_calls && 
		iteration_ctr < 1000)
	{
		previous_rmsd = powell->rmsd_overlay_data->rmsd;

		powell_minimization(powell->rmsd_overlay_data->translate_rotate_array, 7, protein->parameters.max_iterations, &protein->align_rmsd, &rmsd_align_objective_function, powell, start_time);
	 
		if(fabs(previous_rmsd - powell->rmsd_overlay_data->rmsd)<EPS)
			++iteration_ctr;
		else
			iteration_ctr=0; 

	}

	protein->align_rmsd = powell->rmsd_overlay_data->rmsd;

	if(protein->transform_matrix==NULL)
		protein->transform_matrix = (double *)calloc(14,sizeof(double));

	translate_rotate_array_to_transform_matrix(powell->rmsd_overlay_data->translate_rotate_array, protein->transform_matrix);
	
	free_memory(powell->rmsd_overlay_data->initial_coord);
	free_memory(powell->rmsd_overlay_data->target_coord);
	free_memory(powell->rmsd_overlay_data->translate_rotate_array);
	free_memory(powell->rmsd_overlay_data);
	free_memory(powell);

	transform_pdb(protein->Template, protein->transform_matrix);
	pdbATOM_to_mini_pdbATOM(protein->Template,protein->mini_Template);

	pdbATOM_to_CHROMOSOME(protein->Template, &protein->final_chr, protein->resparam);

	dummy_varpos = (VARIABLE_POSITION *)calloc(MAX_RESIDUES, sizeof(VARIABLE_POSITION));
	protein->fixed_atoms = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM));
	
            i=1; j=0;
            while(protein->Template[i].seq_position!=ENDFLAG)
            {
                ++j;
                dummy_varpos[j].seq_position = protein->Template[i].seq_position;
                while(protein->Template[i].seq_position == dummy_varpos[j].seq_position)
                    ++i;

            }
            dummy_varpos[j+1].seq_position = ENDFLAG;

            make_fixed_pdbATOM(protein->Template, dummy_varpos, protein->fixed_atoms);

	protein->mini_fixed_atoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
	pdbATOM_to_mini_pdbATOM(protein->fixed_atoms,protein->mini_fixed_atoms);

	final_chr_to_final_pdb_final_energy(protein);

}