/* EGAD: scmf.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 self-consistent-mean-field optimization (scmf)
   
*/



#include "scmf.h"

/* convert probabs in lookupEnergy.P entries (from scmf optimization) to residue and rotamer probabilities */
void convert_scmf_probabs_to_rotamer_probabs(VARIABLE_POSITION *varPos)
{
    int i, i_res, i_res_rot;
    double rotamer_probab_sum, residue_probab_sum;
        
    i=1;
    while(varPos[i].seq_position!=ENDFLAG)
    {
	residue_probab_sum = 0;
	for(i_res=1;i_res<=varPos[i].number_of_choices;++i_res)
	{
	    rotamer_probab_sum = 0;
	    for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		rotamer_probab_sum += varPos[i].choice[i_res].lookup_res_ptr->lookupRot[i_res_rot].P[0];
	    
	    if(rotamer_probab_sum !=0)
	    {
		for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		    varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].freq = 
			    varPos[i].choice[i_res].lookup_res_ptr->lookupRot[i_res_rot].P[0]/rotamer_probab_sum;
	   
	    
		varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[0] = 0;
		varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[1] = 0;
		varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[2] = 0;
	    
		for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		    varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[i_res_rot] = 
				varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].freq +
				varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[i_res_rot-1];
	     
		residue_probab_sum += varPos[i].lookup_energy_ptr->lookupRes[i_res].P;
	     }
	     
	}
	
	for(i_res=1;i_res<=varPos[i].number_of_choices;++i_res)
	    varPos[i].choice[i_res].composition = varPos[i].lookup_energy_ptr->lookupRes[i_res].P/residue_probab_sum;
	 
	varPos[i].residue_freqs[0] = 0;  varPos[i].residue_freqs[1] = 0; varPos[i].residue_freqs[2] = 0;
	for(i_res=1;i_res<=varPos[i].number_of_choices;++i_res)  
	    varPos[i].residue_freqs[i_res] = varPos[i].choice[i_res].composition + varPos[i].residue_freqs[i_res-1];
	
	++i;
    }
    
}

/* calculates the difference in rotamer probabilties between the current and 
    last SCMF cycle; used for determining scmf convergence */
double conf_matrix_diff(LOOKUP_ENERGY *lookupEnergy, VARIABLE_POSITION *varPos)
{
    int i, i_res, i_res_rot;
    double diff, subdiff;
    
    diff = 0;
    
    i=1;
    while(varPos[i].seq_position!=ENDFLAG)
    {
	for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res)
	{
	    for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
	    {
		subdiff = ERESROT.P[0] - ERESROT.P[1];
		diff += subdiff*subdiff;
	    }
	}
	++i;
    }
    
    return(sqrt(diff));
    
}

/* some structures used locally for scmf */ 
typedef struct
{
  double *rotamer_freq;  
} mini_CHOICE;
 
typedef struct
{
  int seq_position;
  mini_CHOICE *choice;
  double *residue_freq;
} mini_VARIABLE_POSITION; 
 
/* self consistent mean field optimization. 
	assumes that protein has been through input_stuff and generate_lookup_table 
*/

void scmf(PROTEIN *protein)
{
    int i, j, z;
    int i_res, j_res;
    int i_res_rot, j_res_rot;
    double boltz_denominator;
    int seed_flag, convergence_ctr;
    double TdS;
    double convergence_diff;
    double highest_P, lowest_E;
    double E_scmf;
    GENE i_gene, j_gene;
    int best_i_res, best_i_res_rot, ctr;
    int k, numVarPositions, novel_flag;
    CHROMOSOME test_chr;
    extern float NON_INTERACT_PTR;
    extern double RT, MAX_OPTIMIZATION_TIME, SCMF_TEMPERATURE;
    extern LOOKUP_ENERGY_ROTAMER_X *NON_INTERACT_LOOKUP_ROT_X;
    extern LOOKUP_ENERGY_RESIDUE_X *NON_INTERACT_LOOKUP_RES_X;
    double one_minus_lambda, lambda;
    double *past_energies, E_prequench;
    double E_temp, E_pre_quench, rt;
    mini_VARIABLE_POSITION *mini_varpos;
    extern int LOGFILE_FLAG;
    extern char *LOOKUP_TABLE_DIRECTORY;
    FILE *file_ptr;
    time_t now, start_time;
    char *logfile=NULL;
	
	rt = RT;
	RT = SCMF_TEMPERATURE*R_univ;

	strcpy(LOOKUP_TABLE_DIRECTORY,protein->parameters.lookup_energy_table_directory);

	if(LOGFILE_FLAG != 0)
	{
	logfile = (char *)calloc(MAXLINE,sizeof(char));
	sprintf(logfile, "%s.log", protein->parameters.output_prefix);
	file_ptr = fopen_file(logfile, "a");
	fclose(file_ptr);
	}
	
    convergence_ctr = 1; TdS = 0; lowest_E = DBL_MAX; E_scmf = DBL_MAX;
    best_i_res = ENDFLAG; best_i_res_rot = ENDFLAG;

	/* allocate memory */

    protein->E_scmf = DBL_MAX;
    i=1;
    while(protein->var_pos[i].seq_position!=ENDFLAG)
      ++i;
    numVarPositions = i-1;
    
    i_gene = (GENE)malloc(sizeof(MENDEL));
    j_gene = (GENE)malloc(sizeof(MENDEL));
       
    mini_varpos = (mini_VARIABLE_POSITION *)calloc(numVarPositions+2, sizeof(mini_VARIABLE_POSITION));
    for(i=1;i<=numVarPositions;++i)
      {
        mini_varpos[i].seq_position = protein->var_pos[i].seq_position;
	mini_varpos[i].residue_freq = (double *)calloc(protein->var_pos[i].number_of_choices+2, sizeof(double));
	mini_varpos[i].choice = (mini_CHOICE *)calloc(protein->var_pos[i].number_of_choices+2, sizeof(mini_CHOICE));
	for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
	  mini_varpos[i].choice[i_res].rotamer_freq = (double *)calloc(protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers+2, sizeof(double));
      }
    
    
	/* seed the mean-field with probabilities; default to uniform
		weighted uses rotamer probabs; random is completely random */

    seed_flag = 1;
    if(strcmp(protein->parameters.scmf_seed, "unbiased") ==0 || strcmp(protein->parameters.scmf_seed, "UNBIASED") ==0)
      seed_flag = 1;
    else if(strcmp(protein->parameters.scmf_seed, "random") ==0 || strcmp(protein->parameters.scmf_seed, "RANDOM") ==0)
      seed_flag = 2;
    else if(strcmp(protein->parameters.scmf_seed, "weighted") ==0 || strcmp(protein->parameters.scmf_seed, "WEIGHTED") ==0)	
      seed_flag = 3;
    

	start_time = time(NULL);
	z=1;
    while(z<=protein->parameters.number_scmf_cycles)
      {
	for(i=1;i<=numVarPositions;++i)
	  {
	    for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
	      {
		protein->lookupEnergy[i].lookupRes[i_res].P = 0;
		protein->lookupEnergy[i].lookupRes[i_res].TdS = 0;
		protein->lookupEnergy[i].lookupRes[i_res].E_mf = 0;
		
		for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		  {
		    if(protein->ERESROT.P == NULL)
		      protein->ERESROT.P = (double *)calloc(2, sizeof(double));
		    
		    if(seed_flag == 1)
		      protein->ERESROT.P[0] = 1.0;
		    else
		      if(seed_flag == 2)
			protein->ERESROT.P[0] = dice();
		      else
			  {
					protein->ERESROT.P[0] = protein->ERESROT.energy_var_fix;
			
					if(protein->var_pos[i].choice[i_res].in_use_flag==0)	/* penalize not in use choices */
						protein->ERESROT.P[0] += 1e10;
			   }	

		    if(i_res == 1 && i_res_rot == 1)
		      lowest_E = protein->ERESROT.P[0];
		    if(protein->ERESROT.P[0] < lowest_E)
		      lowest_E = protein->ERESROT.P[0];
		  }
	      }
	    
	    boltz_denominator = 0; 
	    for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
	      for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		{
		  protein->ERESROT.P[0] = exp(-(protein->ERESROT.P[0] - lowest_E)/RT);
		  boltz_denominator += protein->ERESROT.P[0]; 
		}
	    
	    for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)  
	      { 
		for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		  {
		    protein->ERESROT.P[0] = protein->ERESROT.P[0]/boltz_denominator;
		  }
	      }
	  }	

	/* parameters for combining current and previous probabilties for generating new probabilties */
	lambda = protein->parameters.scmf_lambda;
	one_minus_lambda =  1.0 - lambda;
		
	/* actual mean field optimization */
	convergence_ctr = 1; convergence_diff = 1;  
	
	while(convergence_diff >= EPS && convergence_ctr <= protein->parameters.max_scmf_iterations)
	  { 
	    E_scmf = protein->parameters.fixedatoms_energy; 
	    TdS=0;
	    
	    for(i=1;i<=numVarPositions;++i)
	      {
			/* calculate mean field energy for each resimer */

		for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
		  {
		    protein->lookupEnergy[i].lookupRes[i_res].E_mf = 0;
		    for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		      {

				protein->ERESROT.P[1] = protein->ERESROT.energy_var_fix;

				if(protein->var_pos[i].choice[i_res].in_use_flag==0)  /* penalize not in use choices */
						protein->ERESROT.P[1] += 1e10;

				protein->lookupEnergy[i].lookupRes[i_res].E_mf += protein->ERESROT.P[0]*protein->ERESROT.energy_var_fix;
			
			for(j=(i+1);j<=numVarPositions;++j)
			  {
			    if(protein->ERESROT.lookupX[j-i].lookupResX != NON_INTERACT_LOOKUP_RES_X)
			      {
				for(j_res=1;j_res <= protein->var_pos[j].number_of_choices;++j_res)
				  {
				    if(protein->ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X)
				      {
						if(protein->ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX == NULL)
						{
							/* make working genes */
      							i_gene->choice_ptr = &protein->var_pos[i].choice[i_res];
      							i_gene->seq_position = protein->var_pos[i].seq_position;
      							i_gene->j_choice_index = i_res-1;
      							i_gene->lookupRot = &i_gene->choice_ptr->lookup_res_ptr->lookupRot[1];
	  						i_gene->lookupRot_index = 1;
							i_gene->varpos_ptr=&protein->var_pos[i];

      							j_gene->choice_ptr = &protein->var_pos[j].choice[j_res];
      							j_gene->seq_position = protein->var_pos[j].seq_position;
      							j_gene->j_choice_index = j_res-1;
							j_gene->lookupRot = &j_gene->choice_ptr->lookup_res_ptr->lookupRot[1];
	      						j_gene->lookupRot_index = 1;
							j_gene->varpos_ptr=&protein->var_pos[j];

      							i_gene->nextgene = NULL; j_gene->nextgene = NULL; 

							load_calc_save_lookupResX(i_gene, j_gene);

						}


					for(j_res_rot=1; j_res_rot <= protein->var_pos[j].choice[j_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++j_res_rot)
					  {
					    if(protein->ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX[j_res_rot-1].energy_var_var != NON_INTERACT_PTR)
					      {
						E_temp = protein->ERESROTj.P[0]*protein->ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX[j_res_rot-1].energy_var_var;
						protein->ERESROT.P[1] += E_temp;
						protein->lookupEnergy[i].lookupRes[i_res].E_mf += protein->ERESROT.P[0]*E_temp;
					      }
					  }
				      }
				  }
			      }
			  }
			
			for(j=(i-1);j>=1;--j)
			  {
			    for(j_res=1;j_res <= protein->var_pos[j].number_of_choices;++j_res)
			      {
				for(j_res_rot=1; j_res_rot <= protein->var_pos[j].choice[j_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++j_res_rot)    
				  {
				    if(protein->ERESROTj.lookupX[i-j].lookupResX != NON_INTERACT_LOOKUP_RES_X)
				      {
					if(protein->ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X)
					  {
						if(protein->ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX == NULL)
						{
							/* make working genes */
      							i_gene->choice_ptr = &protein->var_pos[i].choice[i_res];
      							i_gene->seq_position = protein->var_pos[i].seq_position;
      							i_gene->j_choice_index = i_res-1;
							i_gene->lookupRot = &i_gene->choice_ptr->lookup_res_ptr->lookupRot[1];
	  						i_gene->lookupRot_index = 1;
							i_gene->varpos_ptr=&protein->var_pos[i];
      
      							j_gene->choice_ptr = &protein->var_pos[j].choice[j_res];
      							j_gene->seq_position = protein->var_pos[j].seq_position;
      							j_gene->j_choice_index = j_res-1;
							j_gene->lookupRot = &j_gene->choice_ptr->lookup_res_ptr->lookupRot[1];
	      						j_gene->lookupRot_index = 1;
							j_gene->varpos_ptr=&protein->var_pos[j];

      							i_gene->nextgene = NULL; j_gene->nextgene = NULL; 

							load_calc_save_lookupResX(j_gene, i_gene);

						}



					    if(protein->ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX[i_res_rot-1].energy_var_var != NON_INTERACT_PTR)
					      { 
						E_temp = protein->ERESROTj.P[0]*protein->ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX[i_res_rot-1].energy_var_var;    
						protein->ERESROT.P[1] += E_temp;
						protein->lookupEnergy[i].lookupRes[i_res].E_mf += protein->ERESROT.P[0]*E_temp;
					      }
					  }
				      }
				  }
			      }
			  }
			
			if(i_res == 1 && i_res_rot == 1)
			  lowest_E = protein->ERESROT.P[1];
			if(protein->ERESROT.P[1] < lowest_E)
			  lowest_E = protein->ERESROT.P[1];
			
		      }
		  }
		
			/* update probabilities for each resimer */		

		boltz_denominator = 0;
		for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
		  {
		    for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		      {
			protein->ERESROT.P[1] = exp(-(protein->ERESROT.P[1] - lowest_E)/RT);
			boltz_denominator += protein->ERESROT.P[1];
		      }
		  }
		
		
		for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
		  {
		    protein->lookupEnergy[i].lookupRes[i_res].P = 0;
		    protein->lookupEnergy[i].lookupRes[i_res].TdS = 0;
		    E_scmf += protein->lookupEnergy[i].lookupRes[i_res].E_mf;
		    for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		      {
			/* memory of previous cycle to prevent oscillations per Koehl & Delarue JMB 239: 249 (1994) */
			protein->ERESROT.P[1] = lambda*protein->ERESROT.P[1]/boltz_denominator + one_minus_lambda*protein->ERESROT.P[0];
			protein->lookupEnergy[i].lookupRes[i_res].P += protein->ERESROT.P[1];
			if(protein->ERESROT.P[1] >= EPS)
			  protein->lookupEnergy[i].lookupRes[i_res].TdS += protein->ERESROT.P[1]*log(protein->ERESROT.P[1]);
		      }
		    protein->lookupEnergy[i].lookupRes[i_res].TdS = -RT*protein->lookupEnergy[i].lookupRes[i_res].TdS;
		    TdS += protein->lookupEnergy[i].lookupRes[i_res].TdS;
		  }

	      } /* end i */
	     
	    /* check for convergence */
	    convergence_diff = conf_matrix_diff(protein->lookupEnergy, protein->var_pos);
	    
	
	
	if(LOGFILE_FLAG != 0)
	{         
	file_ptr = fopen_file(logfile, "a");
	now = time(NULL);
	fprintf(file_ptr,"scmf\t%d\t%f\t%f\t%f\t%f\t%s",
					convergence_ctr, convergence_diff, E_scmf, TdS, E_scmf - TdS, ctime(&now) ); 
	fclose(file_ptr);
	}
           
	        
	    
	    for(i=1;i<=numVarPositions;++i)
	      for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
		for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		{
		 	protein->ERESROT.P[0] = protein->ERESROT.P[1];	/* save current probab in P[0] */
	  	}
  
	now = time(NULL);
	if(difftime(now,start_time) >= MAX_OPTIMIZATION_TIME)
	{
		
		file_ptr = fopen_file(logfile, "a");
		now = time(NULL);
		fprintf(file_ptr,"scmf\tout of time at iteration %d\t%s",convergence_ctr, ctime(&now) ); 
		fclose(file_ptr);
		convergence_ctr = protein->parameters.max_scmf_iterations;
	}

	    ++convergence_ctr;

	
	    
	  } /* end mean-field optimization */
	
	/* if doing multiple random seeds, save the lowest-energy scmf soln. */
	if(E_scmf < protein->E_scmf)
	  {
	    
	    protein->E_scmf = E_scmf;
	    
	    protein->final_chr.genes = (MENDEL *)malloc(sizeof(MENDEL));
	    protein->final_chr.firstgene = protein->final_chr.genes;
	    protein->final_chr.genes->nextgene = NULL;
	    protein->final_chr.genes->seq_position = ENDFLAG;
	    	    
		/* generate a structure which has the highest probability rotamer at each position */

	    for(i=1;i<=numVarPositions;++i)
	      {
		mini_varpos[i].residue_freq[0] = 0;
		mini_varpos[i].residue_freq[1] = 0;  
		mini_varpos[i].residue_freq[2] = 0; 
		
		highest_P = 0; 
		for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
		  {
		    
		    mini_varpos[i].residue_freq[i_res] = mini_varpos[i].residue_freq[i_res-1] + protein->lookupEnergy[i].lookupRes[i_res].P;
		    
		    mini_varpos[i].choice[i_res].rotamer_freq[0] = 0; 
		    mini_varpos[i].choice[i_res].rotamer_freq[1] = 0; 
		    mini_varpos[i].choice[i_res].rotamer_freq[2] = 0; 
		    
		    
		    for(i_res_rot=1; i_res_rot<= protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot)
		      {
			if(protein->ERESROT.P[1] > highest_P)
			  {
			    best_i_res = i_res;
			    best_i_res_rot = i_res_rot;
			    highest_P = protein->ERESROT.P[1];
			  }
	
		

			mini_varpos[i].choice[i_res].rotamer_freq[i_res_rot] = mini_varpos[i].choice[i_res].rotamer_freq[i_res_rot-1] + 
			  protein->ERESROT.P[1];
			  
		      }
		     
		  }
		 
		i_res = best_i_res; 
		i_res_rot = best_i_res_rot;
		
		protein->final_chr.genes->seq_position = protein->var_pos[i].seq_position;
		
		protein->final_chr.genes->varpos_ptr = &protein->var_pos[i];
		protein->final_chr.genes->choice_ptr = &protein->var_pos[i].choice[i_res];
		protein->final_chr.genes->lookupRot_index = i_res_rot;
		protein->final_chr.genes->j_choice_index = i_res-1;
		
		protein->final_chr.genes->lookupRot = &protein->ERESROT;
		
		protein->final_chr.genes->chi = (double *)calloc(protein->final_chr.genes->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi+2, sizeof(double));
		
                protein->final_chr.genes->chi[0] = 0;
		k=1;
		while(k<=protein->final_chr.genes->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi)
		{
		  protein->final_chr.genes->chi[k] = protein->final_chr.genes->lookupRot->rotamer->chi[k];
		  ++k;
		}
		protein->final_chr.genes->chi[k]=0;
		
		
		protein->final_chr.genes->nextgene = (MENDEL *)malloc(sizeof(MENDEL));
		protein->final_chr.genes = protein->final_chr.genes->nextgene;
		protein->final_chr.genes->nextgene = NULL;
		protein->final_chr.genes->seq_position = ENDFLAG;
		
	      }
	    protein->final_chr.genes = protein->final_chr.firstgene;
	  }
	if(convergence_ctr >= protein->parameters.max_scmf_iterations)
      {
	if(LOGFILE_FLAG != 0)
	{         
	file_ptr = fopen_file(logfile, "a");
	now = time(NULL);
	fprintf(file_ptr,"self-consistency not achieved before %d cycles\t%s", protein->parameters.max_scmf_iterations, ctime(&now)); 
	fclose(file_ptr);
	}
      }
	seed_flag = 2;
	++z; 
	
	/* out of time; don't try with another random seed */
	now = time(NULL);
	if(difftime(now,start_time) >= MAX_OPTIMIZATION_TIME)
		z = protein->parameters.number_scmf_cycles+10;


      }  /* end all scmf calcs */
    
	RT = rt;

	CHROMOSOME_to_lookupEnergy(&protein->final_chr,&protein->parameters.fixedatoms_energy); 

	/* actually generate a structure */
    if(protein->parameters.scmf_quench_flag == 1)
      {
	/* quench the most-probable-rotamers solution */

	E_prequench = protein->final_chr.energy;
	HQM_rotamers(&protein->parameters, &protein->final_chr, protein->var_pos);
	
	/* if more than one solution is desired */
	if(protein->parameters.number_scmf_solns!=0)
	  {
	    past_energies = (double *)calloc(protein->parameters.number_scmf_solns+2, sizeof(double));
	    for(i=1;i<=protein->parameters.number_scmf_solns;++i)
	      past_energies[i] = 0;
	    
	    past_energies[1] = 	E_prequench;
	    inoculate_sidechains(&test_chr, protein->var_pos,0);

		/* create random structures, biased by the scmf probabilities */
	    for(j=1;j<=protein->parameters.number_scmf_solns;++j)
	      {
			/* generate a unique biased random structure */
		novel_flag=0; ctr=0;
		while(novel_flag == 0 && ctr<=protein->parameters.number_scmf_solns/5)
		  {
		    novel_flag=1;
		    
		    i=1;
		    test_chr.genes = test_chr.firstgene;  
		    while(test_chr.genes->seq_position!=ENDFLAG)    /* re-initialize */ 
		      { 
			/* pick the residuetype */	
			if(protein->var_pos[i].number_of_choices==1)
			  test_chr.genes->j_choice_index = 1; 	   
			else	
			  test_chr.genes->j_choice_index = 
			    roll_loaded_dice(mini_varpos[i].residue_freq, protein->var_pos[i].number_of_choices);
			
			test_chr.genes->choice_ptr = &(protein->var_pos[i].choice[test_chr.genes->j_choice_index]); 
			
			/* assign the rotamer  */	 
			test_chr.genes->lookupRot_index = 
			  roll_loaded_dice(mini_varpos[i].choice[test_chr.genes->j_choice_index].rotamer_freq, test_chr.genes->choice_ptr->resparam_ptr->rotamerlib_ptr->number_of_rotamers);	
			
			/*  copy library rotamer values  */		
			test_chr.genes->chi = 
			  test_chr.genes->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[test_chr.genes->lookupRot_index].chi;	
			
			/* pointer to lookupRotamer  */	  
			test_chr.genes->lookupRot = &(test_chr.genes->choice_ptr->lookup_res_ptr->lookupRot[test_chr.genes->lookupRot_index]);   	 
			test_chr.genes->j_choice_index = test_chr.genes->j_choice_index - 1;   
			
			test_chr.genes = test_chr.genes->nextgene;	
			++i; 
		      }
		    test_chr.genes = test_chr.firstgene;
		    
			/* score the energy of this biased, random structure */
		    CHROMOSOME_to_lookupEnergy(&test_chr,&protein->parameters.fixedatoms_energy); 
		    
		    for(k=1;k<=j;++k)
		      {
			if(test_chr.energy == past_energies[k])
			  novel_flag=0;
		      }
		    
		    ++ctr;
		  }
		past_energies[j]=test_chr.energy;
		
		  /* use this biased random structure as a starting point for MC */
		MC_rotamers(&protein->parameters, &test_chr, protein->var_pos);
		E_pre_quench = test_chr.energy;
		HQM_rotamers(&protein->parameters, &test_chr, protein->var_pos);
		
	if(LOGFILE_FLAG != 0)
	{	 
	file_ptr = fopen_file(logfile, "a");
	now = time(NULL);
	fprintf(file_ptr, "quench\t%d\t%f\t%f\t%s", j, E_pre_quench, test_chr.energy, ctime(&now));
	fclose(file_ptr);
	}	  
		
	    	/* save the best soln to final_chr */
		if(test_chr.energy <= protein->final_chr.energy)
		{
		  copyCHROMOSOME(protein->final_chr, test_chr);
		}
	      }
	    
	    free_memory(past_energies);
	  }
	
	
      }
    
	/* get the final structures and entropy */
	  
	    protein->E_working = protein->final_chr.energy;
	    
	    final_chr_to_final_pdb_final_energy(protein);
	    protein->final_energy->TdS = TdS - protein->final_energy->TdS;  // protein->final_energy->TdS from final_chr_to_final_pdb_final_energy has TdS of unf state
	

	/* free memory */
    for(i=1;i<=numVarPositions;++i)
      {
	for(i_res=1; i_res<=protein->var_pos[i].number_of_choices; ++i_res)
	  free_memory(mini_varpos[i].choice[i_res].rotamer_freq);
	free_memory(mini_varpos[i].residue_freq);
	free_memory(mini_varpos[i].choice);
      }
    free_memory(mini_varpos);
    free_memory(i_gene); free_memory(j_gene);
    
	if(LOGFILE_FLAG != 0)
		free_memory(logfile);
    
}