/* EGAD: HQM_rotamers.h
   
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 is the header for HQM_rotamers.cpp; also contains the inlined
	function find_best_res_res_rot for identifying the best resimer at a given position

*/



#ifndef HQM_rotamers_header_flag
#define HQM_rotamers_header_flag

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <float.h>
#include <limits.h>
#include <math.h>

#include "structure_types.h"
#include "moremath.h"
#include "GA_utilities.h"
#include "CHROMOSOME_to_lookupEnergy.h"
#include "solubility.h"


int has_this_position_been_used(int k, int *already_used_positions);

/* exhaustive downhill line minimization for the sidechains defined by CHROMOSOME chr; used as a quench at the 
    end of GA or MC. chr must be soluble.
*/
void HQM_rotamers(PARAMETERS *parameters, CHROMOSOME *chr, VARIABLE_POSITION *var_pos);

/* exhaustively place all resimers at position i and identify the best_res, best_res_rot (lowest energy) */
inline double find_best_res_res_rot(int i, VARIABLE_POSITION *var_pos, GENE firstgene,
				double fixed_sasa, double fixed_hphob_sasa, double fixed_charge, double fixed_E_transfer,
				int *best_res, int *best_res_rot)
{
	double best_energy;
	int i_res, i_res_rot;
	MENDEL dummygene;
	GENE i_gene;
	double overall_charge, frac_hydrophobic;
	double sasa_total, sasa_hydrophobic, E_transfer_density, E_transfer;
	double sasa_other, sasa_hphob_other, charge_other, energy, E_transfer_other;
	extern int SOLUBILITY_CUTOFF_FLAG;
	extern double FRACTION_HYDROPHOBIC_SASA_CUTOFF, TRANSFER_FREE_ENERGY_DENSITY_CUTOFF;

	if(var_pos[i].fixed_flag!=0)
		failure_report("ERROR position i for HQM_rotamers.h: find_best_res_res_rot must have fixed_flag=0 (not fixed)","abort");


	

	i_gene = &dummygene;

	overall_charge=0; frac_hydrophobic=0;

	i_gene->seq_position = var_pos[i].seq_position;
	i_gene->varpos_ptr = &var_pos[i];

 	best_energy=DBL_MAX;
	 for(i_res=1;i_res<=var_pos[i].number_of_choices;++i_res)
	    {
		if(var_pos[i].choice[i_res].in_use_flag==1)	/* choice must be in use */
		{		

		i_gene->choice_ptr = &var_pos[i].choice[i_res];
		i_gene->j_choice_index = i_res-1; 
		for(i_res_rot=1;i_res_rot<=var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers;++i_res_rot)
		{
			if(i_gene->choice_ptr->lookup_res_ptr->lookupRot[i_res_rot].in_use_flag==1)
			{
			i_gene->lookupRot_index = i_res_rot;
			i_gene->lookupRot = &(i_gene->choice_ptr->lookup_res_ptr->lookupRot[i_gene->lookupRot_index]);
			i_gene->chi = i_gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[i_gene->lookupRot_index].chi;

			energy = lookup_energy_sasa_charge_other(i, i_gene, firstgene, &sasa_other, &sasa_hphob_other, &charge_other, &E_transfer_other) + iROT->energy_var_fix;
			
			if(energy  - best_energy < EPS)
			{
				

				if(SOLUBILITY_CUTOFF_FLAG == 1)
				{
					overall_charge = fixed_charge + charge_other + i_gene->choice_ptr->resparam_ptr->overall_charge_pH_avg;
					sasa_total = fixed_sasa + sasa_other + i_gene->choice_ptr->lookup_res_ptr->sasa_total;
					sasa_hydrophobic = fixed_hphob_sasa + sasa_hphob_other + i_gene->choice_ptr->lookup_res_ptr->sasa_hphob;
					E_transfer = fixed_E_transfer + E_transfer_other + i_gene->choice_ptr->lookup_res_ptr->E_transfer;

					frac_hydrophobic = sasa_hydrophobic/sasa_total;
					E_transfer_density = E_transfer/sasa_total;
					
					if(frac_hydrophobic <= FRACTION_HYDROPHOBIC_SASA_CUTOFF && 
						charge_within_specifications(overall_charge) == 1 &&
						E_transfer_density <= TRANSFER_FREE_ENERGY_DENSITY_CUTOFF)
					{
			    		best_energy = energy;
			    		*best_res = i_res;
			    		*best_res_rot = i_res_rot;
					}
				}
				else
				{
					best_energy = energy;
			    	*best_res = i_res;
			    	*best_res_rot = i_res_rot;
				}
			}
			}

		}
		}  
	    }

	return(best_energy);

}



inline double find_better_res_res_rot(int i, VARIABLE_POSITION *var_pos, GENE firstgene, 
				double fixed_sasa, double fixed_hphob_sasa, double fixed_charge, double fixed_E_transfer,
				int *best_res, int *best_res_rot)
{
	double best_energy;
	int i_res, i_res_rot;
	MENDEL dummygene;
	GENE i_gene;
	double overall_charge, frac_hydrophobic;
	double sasa_total, sasa_hydrophobic, E_transfer_density, E_transfer;
	double sasa_other, sasa_hphob_other, charge_other, energy, E_transfer_other;
	extern int SOLUBILITY_CUTOFF_FLAG;
	extern double FRACTION_HYDROPHOBIC_SASA_CUTOFF, TRANSFER_FREE_ENERGY_DENSITY_CUTOFF;
        double frac;

	if(var_pos[i].fixed_flag!=0)
		failure_report("ERROR position i for HQM_rotamers.h: find_better_res_res_rot must have fixed_flag=0 (not fixed)","abort");

	i_gene = &dummygene;

	i_gene = firstgene;
	while(i_gene->varpos_ptr->varpos_index != i)
		i_gene = i_gene->nextgene;

	frac = 1000.0/((double)i_gene->choice_ptr->resparam_ptr->rotamerlib_ptr->number_of_rotamers);	

	best_energy = lookup_energy_sasa_charge_other(i, i_gene, firstgene, &sasa_other, &sasa_hphob_other, &charge_other, &E_transfer_other) + iROT->energy_var_fix;
	

	i_gene = &dummygene;

	overall_charge=0; frac_hydrophobic=0;

	i_gene->seq_position = var_pos[i].seq_position;
	i_gene->varpos_ptr = &var_pos[i];

 	
	 for(i_res=1;i_res<=var_pos[i].number_of_choices;++i_res)
	    {
		if(var_pos[i].choice[i_res].in_use_flag==1)	/* choice must be in use */
		{		

		i_gene->choice_ptr = &var_pos[i].choice[i_res];
		i_gene->j_choice_index = i_res-1; 
		for(i_res_rot=1;i_res_rot<=var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers;++i_res_rot)
		{
			if(i_gene->choice_ptr->lookup_res_ptr->lookupRot[i_res_rot].in_use_flag==1)
			if(dice()<=frac)
			{
			i_gene->lookupRot_index = i_res_rot;
			i_gene->lookupRot = &(i_gene->choice_ptr->lookup_res_ptr->lookupRot[i_gene->lookupRot_index]);
			i_gene->chi = i_gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[i_gene->lookupRot_index].chi;

			energy = lookup_energy_sasa_charge_other(i, i_gene, firstgene, &sasa_other, &sasa_hphob_other, &charge_other, &E_transfer_other) + iROT->energy_var_fix;
			
			if(energy  - best_energy < EPS)
			{
				

				if(SOLUBILITY_CUTOFF_FLAG == 1)
				{
					overall_charge = fixed_charge + charge_other + i_gene->choice_ptr->resparam_ptr->overall_charge_pH_avg;
					sasa_total = fixed_sasa + sasa_other + i_gene->choice_ptr->lookup_res_ptr->sasa_total;
					sasa_hydrophobic = fixed_hphob_sasa + sasa_hphob_other + i_gene->choice_ptr->lookup_res_ptr->sasa_hphob;
					E_transfer = fixed_E_transfer + E_transfer_other + i_gene->choice_ptr->lookup_res_ptr->E_transfer;

					frac_hydrophobic = sasa_hydrophobic/sasa_total;
					E_transfer_density = E_transfer/sasa_total;
					
					if(frac_hydrophobic <= FRACTION_HYDROPHOBIC_SASA_CUTOFF && 
						charge_within_specifications(overall_charge) == 1 &&
						E_transfer_density <= TRANSFER_FREE_ENERGY_DENSITY_CUTOFF)
					{
			    		best_energy = energy;
			    		*best_res = i_res;
			    		*best_res_rot = i_res_rot;
					}
				}
				else
				{
					best_energy = energy;
			    	*best_res = i_res;
			    	*best_res_rot = i_res_rot;
				}
			}
			}

		}
		}  
	    }

	return(best_energy);

}



#endif