/* EGAD: CHROMOSOME_to_lookupEnergy.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 CHROMOSOME_to_lookupEnergy.cpp.

It also contains the inlined function lookup_energy_sasa_charge_other, which returns the interaction energy of 
	a rotamer with rotamers at all other positions in the structure.
   
*/


#ifndef CHROMOSOME_to_lookupEnergy_header_flag
#define CHROMOSOME_to_lookupEnergy_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 "energy_functions.h"
#include "pairwise_energy_calc.h"
#include "GENES_pdb.h"
#include "output_stuff.h"



void goldstein_singles_elimination(PARAMETERS *parameters, VARIABLE_POSITION *varpos, LOOKUP_ENERGY *lookupEnergy);

/* This function calculates the energy of a CHROMOSOME chr
   by adding up the pre-calculated variable residue/ rotamer energies
   stored in the LOOKUP_ENERGY structure. The genes in the chromosome must be linked to 
	the appropriate pointers. fixed_energy is calculated by  fixedatoms_energy and adjusted by
	generate_lookup_table to include internal energies of non-moving sidechains   */

void CHROMOSOME_to_lookupEnergy(CHROMOSOME *chr, double *fixed_energy);


/* protein->final_chr contains the final solution CHROMOSOME. This function allocates
	memory for protein->final_energy, protein->final_pdb, and protein->final_minipdb, 
	calculates the coordinates, structures, and energy components. Should be called 
	before output_PROTEIN */
void final_chr_to_final_pdb_final_energy(PROTEIN *protein);


/* a_gene and b_gene are GENEs with all the pointers and values assigned (by mutate_sidechain, inoculate_sidechains)
	this function returns in the pair energy between the two rotamers (NON_INTERACT_PTR if they do not interact).
   If the lookup table for these residues at these positions have not been calculated or read from disk, this function
	takes care of that.
   If the entire table has been precalculated (default for all except for rotamer_foremen jobs), this function is not
	called.
*/
double load_calc_save_lookupResX(GENE a_gene, GENE b_gene);


inline int find_j_res_rot_minus1_in_hbond_rot_pair(int j_res_rot_minus1, LOOKUP_ROTAMER_X_HBOND *lookupRotX_hbond)
{
	int h;
	int lower,upper;


  	lower = 0; upper = lookupRotX_hbond->num_hbond_rot_pairs-1;

	while(lower<=upper)
    	{
      		h = (lower+upper)/2;
      		if(lookupRotX_hbond->hbond_rot_pair[h].j_index == j_res_rot_minus1)
				return(h);  
      		else
        		if(lookupRotX_hbond->hbond_rot_pair[h].j_index < j_res_rot_minus1)
          			lower = h+1;
        		else
          			upper = h-1;
    	}

	return(ENDFLAG);

}

/* returns the interaction energy of the rotamer encoded by i_gene with all other rotamers in the current structure.
	i = var_pos index of i_gene
	firstgene = firstgene of the chromosome containing the current structure
	these are calculated if SOLUBILITY_CUTOFF_FLAG=1 (solubility as design restraint):
		sasa_other = surface area of all other moving non-i residues
		sasa_hphob_other = hydrophobic surface area of all other moving non-i residues
		charge_other = total charge of all other moving non-i residues
		E_transfer_other = water-octanol transfer energy of all other moving non-i residues
*/

inline double lookup_energy_sasa_charge_other(int i, GENE i_gene, GENE firstgene, 
				double *sasa_other, double *sasa_hphob_other, double *charge_other, double *E_transfer_other)
{
	double current_energy;
	int j,i_minus_j, j_minus_i;
	extern float FIXED_POSITION_PTR, NON_INTERACT_PTR;
	extern LOOKUP_ENERGY_ROTAMER_X *NON_INTERACT_LOOKUP_ROT_X;
        extern LOOKUP_ENERGY_RESIDUE_X *NON_INTERACT_LOOKUP_RES_X;
	extern int SOLUBILITY_CUTOFF_FLAG;
	GENE j_gene;
        extern bool *FIXED_HBOND_SATISF_ARRAY;
        extern ATOMRESPARAM **FIXED_ATOM_HBOND_ATOMRESPARAM;
  	extern int FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH, HBOND_SPECIFICITY_FLAG, BKBN_HBOND_STATUS_ARRAYLENGTH;
	int ii,jj,h,ii2,jj2,hbond_status_arraylength,i_res_rot_minus1, j_res_rot_minus1, x,k;
  	static bool *hbond_status_array=NULL;

	if(i_gene->lookupRot->energy_var_fix==FIXED_POSITION_PTR)
		failure_report("ERROR i_gene sent to lookup_energy_sasa_charge_other must be a bona-fide floating position","abort");

	ii=1; jj=1; ii2=1; jj2=1;

	if(SOLUBILITY_CUTOFF_FLAG == 1)
	{
		*sasa_other=0;
		*sasa_hphob_other=0;
		*charge_other=0;
		*E_transfer_other=0;
	}

	hbond_status_arraylength=FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH;

	if(HBOND_SPECIFICITY_FLAG == 1)
	{
		if(hbond_status_array == NULL)
			hbond_status_array = (bool *)calloc(MAX_ATOMS,sizeof(bool));	

		falsify_bool_array(hbond_status_array,1,MAX_ATOMS-1);		

		copy_bool_array(hbond_status_array, 1, 
								FIXED_HBOND_SATISF_ARRAY, 1, FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH);
		
		
		ii = FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH + 1;

			j_gene = firstgene;	
			while(j_gene->seq_position != i_gene->seq_position)
			{
				if(j_gene->lookupRot->energy_var_fix!=FIXED_POSITION_PTR) 
				{
					for(x=0;x<jROT->num_var_fix_hbond;++x)
							hbond_status_array[jROT->bkbn_hbond_indicies[x]] = 1;

					
					if(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
					{
						copy_bool_array(hbond_status_array, ii, 
								jROT->sidechain_hbond_status, 1, j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms);

						ii = ii + j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms;
						hbond_status_arraylength += j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms;
					}
				}

				j_gene = j_gene->nextgene;
			}
		
			if(i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
					{
						for(x=0;x<iROT->num_var_fix_hbond;++x)
							hbond_status_array[iROT->bkbn_hbond_indicies[x]] = 1;


						copy_bool_array(hbond_status_array, ii, 
								iROT->sidechain_hbond_status, 1, i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms);

						hbond_status_arraylength += i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms;
					}
			ii2 = ii+i_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1;


			j_gene = j_gene->nextgene;

			jj = ii + i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms;
			while(j_gene->seq_position != ENDFLAG)
			{
				if(j_gene->lookupRot->energy_var_fix!=FIXED_POSITION_PTR)
				{
					for(x=0;x<jROT->num_var_fix_hbond;++x)
							hbond_status_array[jROT->bkbn_hbond_indicies[x]] = 1;

		
					if(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
					{
						copy_bool_array(hbond_status_array, jj, 
								jROT->sidechain_hbond_status, 1, j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms);

						jj = jj + j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms;
						hbond_status_arraylength += j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms;
					}

				}
				j_gene = j_gene->nextgene;
			}

			jj = FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH + 1;
		
	}

	
	current_energy = 0;	

        if(i_gene->choice_ptr->in_use_flag==0) /* penalize choices not in use */
		current_energy += 1e10;

	j_gene = firstgene; j=1;
	while(j_gene->seq_position!=ENDFLAG)
	{
	    if(j_gene->lookupRot->energy_var_fix!=FIXED_POSITION_PTR) 
	    {	
		if(HBOND_SPECIFICITY_FLAG == 1)
		{
			if(j_gene->seq_position!=i_gene->seq_position)
			{
				if(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
					jj2 = jj+j_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1;
			}
			else
				if(i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
					jj2 = jj+i_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1;
		}
		
	    if(j_gene->seq_position > i_gene->seq_position)
	    {
		j_minus_i = j-i;

		if(i_gene->lookupRot->lookupX[j_minus_i].lookupResX != NON_INTERACT_LOOKUP_RES_X)
		    if(i_gene->lookupRot->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX!= NON_INTERACT_LOOKUP_ROT_X)
			{
				if(i_gene->lookupRot->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX == NULL)
					load_calc_save_lookupResX(i_gene, j_gene);

				j_res_rot_minus1 = j_gene->lookupRot_index - 1;
				 
					if(i_gene->lookupRot->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX[j_res_rot_minus1].energy_var_var != NON_INTERACT_PTR)
					{
			    		current_energy += 
	     i_gene->lookupRot->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX[j_res_rot_minus1].energy_var_var;

						if(HBOND_SPECIFICITY_FLAG == 1)
						{
							{


							if(i_gene->lookupRot->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX_hbond!=NULL)
							{
							
								h=find_j_res_rot_minus1_in_hbond_rot_pair(j_res_rot_minus1, i_gene->lookupRot->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX_hbond);
			  		   
 								if(h!=ENDFLAG)
								{
									union_bool_arrays(hbond_status_array, ii, hbond_status_array, ii,  ii2,
	iROT->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX_hbond->hbond_rot_pair[h].i_hbond_status, 
								0, i_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1);

									union_bool_arrays(hbond_status_array, jj, hbond_status_array, jj,  jj2,
	iROT->lookupX[j_minus_i].lookupResX[j_gene->j_choice_index].lookupRotX_hbond->hbond_rot_pair[h].j_hbond_status, 
								0, j_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1);

								}
							}
							}
						}

					}
			} 
	    }
	    else if(j_gene->seq_position < i_gene->seq_position)
	    {
			i_minus_j = i-j;

		if(j_gene->lookupRot->lookupX[i_minus_j].lookupResX != NON_INTERACT_LOOKUP_RES_X)
		    if(j_gene->lookupRot->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX != NON_INTERACT_LOOKUP_ROT_X)
			{
				if(j_gene->lookupRot->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX == NULL)
					load_calc_save_lookupResX(j_gene, i_gene);

			i_res_rot_minus1 = i_gene->lookupRot_index-1;


					if(j_gene->lookupRot->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX[i_res_rot_minus1].energy_var_var != NON_INTERACT_PTR)
					{
			    		current_energy += 
		j_gene->lookupRot->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX[i_res_rot_minus1].energy_var_var;


						if(HBOND_SPECIFICITY_FLAG == 1)
						{
							{
					
						if(j_gene->lookupRot->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX_hbond!=NULL)
							{

								h = find_j_res_rot_minus1_in_hbond_rot_pair(i_res_rot_minus1, j_gene->lookupRot->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX_hbond);
			  		   
    								if(h!=ENDFLAG)
								{
									union_bool_arrays(hbond_status_array, ii, hbond_status_array, ii,  ii2,
	jROT->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX_hbond->hbond_rot_pair[h].j_hbond_status, 
								0, i_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1);

									union_bool_arrays(hbond_status_array, jj, hbond_status_array, jj,  jj2,
	jROT->lookupX[i_minus_j].lookupResX[i_gene->j_choice_index].lookupRotX_hbond->hbond_rot_pair[h].i_hbond_status, 
								0, j_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1);

								}

							}
							}
						}


					}
			}
	    }
    
	    if(SOLUBILITY_CUTOFF_FLAG == 1)
	    {
	    	if(j_gene->seq_position != i_gene->seq_position)
				{
	    				*sasa_other += j_gene->choice_ptr->lookup_res_ptr->sasa_total;
	    				*sasa_hphob_other += j_gene->choice_ptr->lookup_res_ptr->sasa_hphob;
					*charge_other += j_gene->choice_ptr->resparam_ptr->overall_charge_pH_avg;
					*E_transfer_other += j_gene->choice_ptr->lookup_res_ptr->E_transfer;
				}	
	    }

		if(HBOND_SPECIFICITY_FLAG == 1)
		{
			if(j_gene->seq_position!=i_gene->seq_position)
			{
				if(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
					jj = jj2 + 1;
			}
			else 
				if(i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
				{
					jj = jj2 + 1;
				}
		}

	    }
	    j_gene = j_gene->nextgene; 
	    ++j;
	}

	if(HBOND_SPECIFICITY_FLAG == 1)
	{
	    for(h=1;h<=BKBN_HBOND_STATUS_ARRAYLENGTH ;++h)
		current_energy -= hbonding_status_dependent_specificity_energy(hbond_status_array[h],FIXED_ATOM_HBOND_ATOMRESPARAM[h]);

	    ii = BKBN_HBOND_STATUS_ARRAYLENGTH+1;
	    j_gene = firstgene;
	    while(j_gene->seq_position!=ENDFLAG)
	    {
		if(j_gene->lookupRot->energy_var_fix==FIXED_POSITION_PTR)
		{
		if(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
		{
			ii2 = ii + j_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1;


			if(any_false(hbond_status_array,ii,ii2)==1)
				current_energy -= j_gene->choice_ptr->resparam_ptr->E_specificity_core;

			if(j_gene->choice_ptr->resparam_ptr->ligand_flag==1)
			{
				k=1;
				for(h=ii;h<=ii2;++h)
				{
					current_energy -= hbonding_status_dependent_specificity_energy(hbond_status_array[h],j_gene->choice_ptr->resparam_ptr->sidechain_hbond_atom_ptr[k]);
					++k;
				}
			}

			ii = ii2+1;
		}
		}
		j_gene = j_gene->nextgene; 
	    }
	 	


	    ii=FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH+1;
	    j_gene = firstgene;
	    while(j_gene->seq_position!=ENDFLAG)
	    {
		if(j_gene->lookupRot->energy_var_fix!=FIXED_POSITION_PTR)
		{
		if(j_gene->seq_position!=i_gene->seq_position)
		{
			if(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
			{

				ii2 = ii + j_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1;

				if(any_false(hbond_status_array,ii,ii2)==1)
					current_energy -= j_gene->choice_ptr->resparam_ptr->E_specificity_core;

				
				if(j_gene->choice_ptr->resparam_ptr->ligand_flag==1)
				{
				k=1;
				for(h=ii;h<=ii2;++h)
				{
					current_energy -= hbonding_status_dependent_specificity_energy(hbond_status_array[h],j_gene->choice_ptr->resparam_ptr->sidechain_hbond_atom_ptr[k]);
					++k;
				}
				}

				ii = ii2+1;

			}

		}
		else
		{
			if(i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms!=0)
			{

				ii2 = ii + i_gene->choice_ptr->resparam_ptr->numhbondsideatoms_minus1;

				if(any_false(hbond_status_array,ii,ii2)==1)
					current_energy -= i_gene->choice_ptr->resparam_ptr->E_specificity_core;
				
				if(i_gene->choice_ptr->resparam_ptr->ligand_flag==1)
				{
				k=1;
				for(h=ii;h<=ii2;++h)
				{
					current_energy -= hbonding_status_dependent_specificity_energy(hbond_status_array[h],i_gene->choice_ptr->resparam_ptr->sidechain_hbond_atom_ptr[k]);
					++k; 
				}
				}

				ii = ii2+1;

			}	
		}
		}	
		j_gene = j_gene->nextgene; 
	    }

	
	    

	/*
	if(DEBUG == 2)
		{
 	printf("partial_structure: %d\t%d\t%d\n",count_false(hbond_status_array,1,hbond_status_arraylength),hbond_status_arraylength,FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH);
	 	fprintf_bool_array(stdout, hbond_status_array,1,hbond_status_arraylength);
		}
	*/	

	}

	return(current_energy);
}


#endif