/* EGAD: VARIABLE_POSITION.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 setting up VARIABLE_POSITION and parsing the user-defined inputfiles.

These functions are called by modify_VARIABLE_POSITION and by input_stuff.cpp: input_VARIABLE_POSITIONS

*/

#include "VARIABLE_POSITION.h"

/* Given an array var_seqpos of variable positions (ending with ENDFLAG), a pdb structure pdb, a distance cutoff, 
   and a scale factor window, this function returns the number of neighbors and variable positions; 
   these are also included in the array neighbors.    
*/
int find_neighbors(int *var_seqpos, pdbATOM *pdb, int *neighbors, double *cutoff, double window)
{
  pdbATOM *residueXatoms;
  int i,j,k,n,flag;
  int pdb_ctr, resXatoms_ctr;
  int numberofNeighbors;
  int *tempout, *tempout2;
  int dummyint, temp_seq;
  int first,last;
  double real_cutoff;
  
  first=1;

  tempout = (int *)calloc(MAX_ATOMS*10+10, sizeof(int));
  tempout2  = (int *)calloc(MAX_ATOMS*10+10, sizeof(int));

  for(i=0;i< MAX_ATOMS*10+10;++i)
  {
	tempout[i] = 0;
  	tempout2[i] = 0;
  }


  residueXatoms = (pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM));	
  temp_seq = var_seqpos[1];
    


  i=1; j=1; 
  while(var_seqpos[i]!=ENDFLAG)
    {
	k=1;
      /* generate residueXatoms for var_seqpos[i] */
      pdb_ctr=1; resXatoms_ctr=1;
      while(pdb[pdb_ctr].seq_position!=var_seqpos[i] && pdb[pdb_ctr].atom_number!=ENDFLAG)  
	++pdb_ctr;    /* advance past preceeding residues */
      while(pdb[pdb_ctr].seq_position==var_seqpos[i] && pdb[pdb_ctr].atom_ptr->other_info<0 && 
	    pdb[pdb_ctr].atom_ptr->other_info!=8 && pdb[pdb_ctr].atom_number!=ENDFLAG)
	++pdb_ctr;   /* advance past backbone for residue var_seqpos[i] */
      while(pdb[pdb_ctr].seq_position==var_seqpos[i] && pdb[pdb_ctr].atom_number!=ENDFLAG)
	{  /* copy sidechain atoms for residue number var_seqpos[i] */
	  residueXatoms[resXatoms_ctr] = pdb[pdb_ctr];  
	  ++resXatoms_ctr; ++pdb_ctr;
	}
      strcpy(residueXatoms[resXatoms_ctr].residuetype,"END");
      resXatoms_ctr=1;
      while(strcmp(residueXatoms[resXatoms_ctr].residuetype,"END")!=0)
	{
	  pdb_ctr=1;
	  while(pdb[pdb_ctr].atom_number!=ENDFLAG)
	    {
	      if(pdb[pdb_ctr].atom_ptr->other_info>0 && strcmp(pdb[pdb_ctr].residuetype,"GLY")!=0 && 
		 strcmp(pdb[pdb_ctr].residuetype,"PRO")!=0 && strcmp(pdb[pdb_ctr].residuetype,"CYD")!=0)
		/*  don't use the backbone atoms, and since gly, pro,disulfide are always fixed, don't count them */
		{
		  dummyint = pdb[pdb_ctr].seq_position;
		  if((*cutoff)>0)	/* use the inputed cutoff distance */
		    real_cutoff=(*cutoff)*(*cutoff);
		  else
		    if((*cutoff)==0)	/* vdw contact */
		      real_cutoff = 
			(double)pow(window*(residueXatoms[resXatoms_ctr].atom_ptr->atom_ptr->sigma/2.0 + pdb[pdb_ctr].atom_ptr->atom_ptr->sigma/2.0), 2.0);
		    else
		      real_cutoff=0;
		  if(Distance_sqrd(residueXatoms[resXatoms_ctr].coord, pdb[pdb_ctr].coord)<=real_cutoff)
		    {
					
		      while(pdb[pdb_ctr].seq_position==dummyint && pdb[pdb_ctr].atom_ptr!=NULL) 
			++pdb_ctr; /* move past this residue since we've already marked it */
		      
		      if(strcmp(pdb[pdb_ctr-1].residuetype,"GLY")!=0 && strcmp(pdb[pdb_ctr-1].residuetype,"NTE")!=0 &&
			 strcmp(pdb[pdb_ctr-1].residuetype,"PRO")!=0 && strcmp(pdb[pdb_ctr-1].residuetype,"CTE")!=0 && 
			 strcmp(pdb[pdb_ctr-1].residuetype,"CYD")!=0)   
			{	
				tempout2[k]=pdb[pdb_ctr-1].seq_position;
			    ++k;

			}
		    }
		  else 
		    ++pdb_ctr;
		}
	      else
		++pdb_ctr; 
	    }
	  ++resXatoms_ctr;
	}

	
	/* sort and compress the list of neighboring residues */
  last=k-1;
  sort_int(&first,&last,tempout2);
  tempout2[k]=ENDFLAG;  /* mark the end of the array */
  
  k=1; 
  while(tempout2[k]!=ENDFLAG)
    {
	flag=0;
      for(n=1;n<j;++n)
	{
	   if(tempout[n]==tempout2[k])	/* we've already got this residue */
		flag = 1;
	   while(tempout[n]==tempout2[k]) /* scroll past the same resnums in the list */
		++k;	
	}	

      if(flag == 0)
      {	      
      	tempout[j]=tempout2[k];
	while(tempout[j]==tempout2[k]) /* scroll past the same resnums in the list */
		++k;
      	++j;
      }
	
    }

	
      ++i;
    }
  
  last=j-1;
  
  /* sort and compress the list of neighboring residues */
  sort_int(&first,&last,tempout);
  tempout[j]=ENDFLAG;  /* mark the end of the array */
  
  i=1;j=1;
  while(tempout[j]!=ENDFLAG)
    {
      neighbors[i]=tempout[j];
      while(neighbors[i]==tempout[j]) /* scroll past the same resnums in the list */
	++j;
      ++i;
    }
  neighbors[i]=ENDFLAG; /* mark the end of the array */
  numberofNeighbors = i-1;
  
  free_memory(tempout); free_memory(tempout2); 
  free_memory(residueXatoms);	  

  return numberofNeighbors;
}

/* This function binds a group of VARIABLE_POSITION arrays. 
    fixed_res contains the positions that are fixed in rotamer and sequence, but are allowed to undergo 
	local rotamer relaxation.
    neighbors contains the positions that are fixed in sequence,  but are allowed to sample all rotamers. 
    user_defined contains the residues that are defined by the user as being variable.
    The array melded is the final VARIABLE_POSITION. 
*/
void meld_VARIABLE_POSITION(VARIABLE_POSITION *melded, VARIABLE_POSITION *fixed_res, 
			    VARIABLE_POSITION *neighbors, VARIABLE_POSITION *user_defined, int *fixed_positions)
{
  int i,j, k, n, flag2, first, numrot, flag, flag3;
  BACKBONE bkbn;
  ROTAMERLIB dummy_rotamerlib;
  double denom, max_freq;
  
  first=1;
  
  i=1;
  while(fixed_res[i].seq_position != ENDFLAG)  /* melded varpos array has the native residue and conf unless specified otherwise */
    {
      melded[i] = fixed_res[i];
      melded[i].fixed_flag=1;
      
      ++i;
    }
  melded[i].seq_position = ENDFLAG;
  
  
  i=1;  /* copy user_defined varpos to over-ride the the options for neighbors */
  while(neighbors[i].seq_position != ENDFLAG)
    {
      j=1; flag2 = 0;
      while(user_defined[j].seq_position != ENDFLAG && flag2==0 )
	{
          
	  if(user_defined[j].seq_position == neighbors[i].seq_position)   
	    {
	      neighbors[i] = user_defined[j];
	      neighbors[i].fixed_flag=0;
	      flag2=1;
	    }
	  ++j;
	} 
      ++i;
    }  

 i=1; /* for each position, see if we need to over-write the native residue with that stored in neighbors (actual neighbors + user-defined) */ 
  while(melded[i].seq_position != ENDFLAG)
    {
      flag3=0;
      if(fixed_positions != NULL)
      {	
      n=1;
      while(fixed_positions[n]!=ENDFLAG)   /* this position is user-defined fixed, so use only the endogeneous rotamer */
	{
	  if(melded[i].seq_position == fixed_positions[n])
	    flag3 = 1;
	  ++n;
	}
      }
	  
      if(flag3 == 0)
	{
      j=1; flag2=0;
      bkbn = melded[i].bkbn;
      while(neighbors[j].seq_position != ENDFLAG && flag2==0)
	{ 
	    if(neighbors[j].seq_position == melded[i].seq_position )
	    {
	      k=1; flag=0;
	      while(flag==0)	/* look for the natural residuetype in the choices available in neighbors */
		{
		  if(neighbors[j].choice[k].resparam_ptr->residuecode == melded[i].choice[1].resparam_ptr->residuecode)
		    flag=1;
		  else
		    {	
		      ++k;
		      if(k> neighbors[j].number_of_choices)  /* the natural residue is not an option at this position */
			flag=2;    	
		    }	
		}
	      if(flag==1)    /* add natural rotamer if the natural residue is an option for this position  */
		{   
		  dummy_rotamerlib = *neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr;    
		  neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr = (ROTAMERLIB *)malloc(sizeof(ROTAMERLIB));
		  *neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr = dummy_rotamerlib;

			if(neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers>1)
		    {
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers = 
			neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers + 1;

		      numrot = neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers;	    
	
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->freq_array = (double *)calloc((numrot + 3), sizeof(double));
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer = (ROTAMER *)calloc((numrot + 3), sizeof(ROTAMER));

		      max_freq=0; denom=0; 
		      for(n=1;n<=(numrot-1);++n)
			{
			  neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n] = dummy_rotamerlib.rotamer[n];
			  
			  if(dummy_rotamerlib.rotamer[n].freq > max_freq)
			    max_freq = dummy_rotamerlib.rotamer[n].freq;
			  denom += dummy_rotamerlib.rotamer[n].freq;	
			}
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot] = melded[i].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1];
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].freq = max_freq+EPS;
		      denom += max_freq;

		      for(n=1;n<=numrot;++n)
			neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].freq = neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].freq/denom;    

			if(neighbors[j].choice[k].resparam_ptr->ligand_flag==0)
		     sort_ROTAMER(&first,&numrot,neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer); 
			else
			{
				if(neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->numChi==1)
				{
					for(n=1;n<=numrot;++n)
						neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].chi[1] = n;
				}
			}

		      for(n=1;n<=numrot;++n)
			{
			  neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].index = n;
		    
			}
		 
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->freq_array[0] = 0; 
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->freq_array[1] = 0;
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->freq_array[2] = 0;
		      for(n=1;n<=numrot;++n)
			{
			  neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].rank = n;
			  neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->freq_array[n] = 
			    neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].freq + 
			    neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->freq_array[n-1];
			}
		      neighbors[j].choice[k].resparam_ptr->rotamerlib_ptr->freq_array[numrot+1] = ENDFLAG;
		    }
		}
	      free_memory(melded[i].choice[1].resparam_ptr);	      
	      melded[i] = neighbors[j];
	      melded[i].bkbn = bkbn;
	      flag2=1;
	    }
	  ++j;
	}
	}
	
      ++i;   
    }
  
    	

}



/* This function generates a VARIABLE_POSITION array for an array listofvarpos of numVarPos variable positions 
	(ending with ENDFLAG) for a structure pdb. 
*/
void listofvarpos_to_VARIABLE_POSITION(int *listofvarpos, VARIABLE_POSITION *var_pos, pdbATOM *pdb, int numVarPos, RESPARAM *resparam)
{
  int i,j, q;
  
  var_pos[numVarPos+1].seq_position = ENDFLAG;
  
  j=1; q=1;
  for(i=1;i<=numVarPos;++i)
    {
      while(pdb[j].seq_position!=listofvarpos[i] && pdb[j].seq_position!=ENDFLAG)   /* advance pdb to next residue */
	++j;
      
      var_pos[i].seq_position = listofvarpos[i];  /* copy seq_position */
      
      while(pdb[j].seq_position!=ENDFLAG && pdb[j].atom_ptr->atomnumber != 8)  /* move to CB */
	++j;
      
      if(pdb[j].seq_position!=ENDFLAG)
	{
	  /* get resparam  pointer assigned */
	  var_pos[i].choice[1].resparam_ptr = (RESPARAM *)malloc(sizeof(RESPARAM));
	  *var_pos[i].choice[1].resparam_ptr = resparam[search_residuetype_RESPARAM(pdb[j].residuetype,resparam)];	
      
	  var_pos[i].number_of_choices = 1;
	  var_pos[i].choice[1].composition = 1.0;
	  var_pos[i].residue_freqs[1] = var_pos[i].choice[1].composition;
	  var_pos[i].residue_freqs[2]=ENDFLAG;
	}
    }
}


/* Given an array varpos containing desired seq positions (ending w/ endflag), this function 
	extracts sidechain dihedrals from pdbATOM array pdb and creates a rotamer library (w/ rotamerlets) 
	for the native conformation for each residue specified in varpos.
	respram = main resparam structure 
*/
void pdbATOM_to_VARIABLE_POSITION(VARIABLE_POSITION *varpos, pdbATOM *pdb, RESPARAM *resparam)
{
  int varpos_index, pdb_index, respdb_index, bkbn_index;
  pdbATOM *residuepdb;
  int q, m, d, n, nn, dd;
  BACKBONE *bkbn;
  ROTAMERLIB rotamerlib;
  
  bkbn = (BACKBONE *)calloc(MAX_RESIDUES, sizeof(BACKBONE));
  residuepdb = (pdbATOM *)calloc((MAX_RES_SIZE + 2), sizeof(pdbATOM));  
  
  extract_bkbn(pdb, bkbn);
  
  varpos_index = 1; pdb_index = 1; bkbn_index = 1;
  while(varpos[varpos_index].seq_position!=ENDFLAG)
    {
      /* advance pdb_index to this variable residue  */
      while(pdb[pdb_index].seq_position!=varpos[varpos_index].seq_position && pdb[pdb_index].seq_position!=ENDFLAG)
	++pdb_index;
      
      while(varpos[varpos_index].seq_position!=bkbn[bkbn_index].seq_position && bkbn[bkbn_index].seq_position!=ENDFLAG)
	++bkbn_index;
      varpos[varpos_index].bkbn = bkbn[bkbn_index];
      
      varpos[varpos_index].choice[1].bkbn = &varpos[varpos_index].bkbn;
      
      varpos[varpos_index].number_of_choices = 1;
      
      /* copy the pdb elements for this variable residue to residuepdb */
      respdb_index = 1;
      while(varpos[varpos_index].seq_position == pdb[pdb_index].seq_position && pdb[pdb_index].seq_position!= ENDFLAG)  
	{
	    
	  residuepdb[respdb_index] = pdb[pdb_index];
	  if(residuepdb[respdb_index].atom_ptr->atomnumber == 8)
	    {
	      varpos[varpos_index].choice[1].resparam_ptr = (RESPARAM *)malloc(sizeof(RESPARAM));
	      
	      *varpos[varpos_index].choice[1].resparam_ptr = resparam[search_residuetype_RESPARAM(residuepdb[respdb_index].residuetype, resparam)];
	      
	    }
	  
	  ++respdb_index; ++pdb_index;
	  
	}
      strcpy(residuepdb[respdb_index].residuetype,"END");
      residuepdb[respdb_index].sasa = ENDFLAG;
      residuepdb[respdb_index].seq_position = ENDFLAG;
      residuepdb[respdb_index].atom_ptr = NULL;
      


      rotamerlib = *varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr;
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr = (ROTAMERLIB *)malloc(sizeof(ROTAMERLIB));
      *varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr = rotamerlib;
      
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->number_of_rotamers = 1;
      
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->freq_array = (double *)calloc(3, sizeof(double));
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->freq_array[0] = 0;
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->freq_array[1] = 1.0;
      
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer = (ROTAMER *)calloc(3, sizeof(ROTAMER));
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].freq = 1.0;
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rank = 1;
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].index = 1;
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].native_rotamer_flag = 1;

	

      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].chi = (double *)calloc(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi+2,sizeof(double));
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].std_dev = (double *)calloc(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi+2,sizeof(double));
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].half_std_dev = (double *)calloc(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi+2,sizeof(double));

      extract_dihedrals(residuepdb, varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].chi, varpos[varpos_index].choice[1].resparam_ptr);

	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].std_dev[0]=0;
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].half_std_dev[0]=0;
      q=1;
      while(q<=varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi)
      {
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].std_dev[q] = NATURAL_ROTAMER_SPREAD;
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].half_std_dev[q] = 0.49*NATURAL_ROTAMER_SPREAD;
        ++q;
      }
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].std_dev[q] = 0;
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].half_std_dev[q] = 0;
      
      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets = 0;
      
      if(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi >= 2)
	{
	  if(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi >= 3)
	    varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets = 27;
	  else    
	    varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets = 9;
	}  	
      else
	varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets = 3;
      
	if(rotamerlib.number_of_rotamers==1)
		varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets = 1;
      
      if(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets != 0)
	{	
	  varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet = (double **)calloc(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets+2, sizeof(double *));
	  for(q=1;q<=varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets;++q)
            {
	      varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q] = (double *)calloc(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi+2, sizeof(double));

	      for(n=0;n<=varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->numChi+1;++n)
		varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][n] = varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].chi[n];

            }
	  
	  
	  if(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets >= 3)
	    {
	      q=1;
	      for(m=-1;m<=1;++m)
		{
		  varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][1] = varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].chi[1] +   
		    m*NATURAL_ROTAMER_SPREAD;
		  if(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets>=9)
		    {
		      n=q;
		      for(d=-1;d<=1;++d)
			{
			  varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][1] = varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[n][1];
			  varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][2] = varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][2]  +  
			    d*NATURAL_ROTAMER_SPREAD;
			  
			  if(varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->num_rotamerlets>=27)
			    {
			      nn = q;
			      for(dd = -1;dd<=1;++dd)
				{
				  varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][1] = varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[nn][1];
				  varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][2] = varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[nn][2];
				  varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][3] = varpos[varpos_index].choice[1].resparam_ptr->rotamerlib_ptr->rotamer[1].rotamerlet[q][3] +
				    dd*NATURAL_ROTAMER_SPREAD;
				  ++q;						    
				}
			    }
			  else	 					    
			    ++q;						    
			}		  		     
		    }
		  else
		    ++q;
		}
	    }
	}
      ++varpos_index;
    }
  
  free_memory(bkbn);  
  free_memory(residuepdb);
}

/* Given a PROTEIN containing PARAMETERS and an array of VARIABLE_POSITION, this 
	function will modify the PARAMETERS and VARIABLE_POSITION inside PROTEIN to include 
	neighboring residues whose rotamer identities are allowed to float, extract the native rotamers
	if applicable, and attach all the appropriate pointers to VARIABLE_POSITIONS. 
	The born radii and energies of fixed atoms are calculated.
	Called by input_stuff.cpp: input_VARIABLE_POSITIONS
*/
void modify_VARIABLE_POSITION(PROTEIN *protein, int *fixed_positions)
{
  int *tempseqVarpos1, *tempseqVarpos, *seqVarpos;
  double max_freq, denom, p;
  int num_residues, original_num_varpos;
  VARIABLE_POSITION *user_defined_varpos, *neighbor_varpos, *fixed_res_VARIABLE_POSITION;
  RESPARAM *fixed_res_RESPARAM;
  CHOICE *neutral_choice;
  ROTAMERLIB dummy_rotamerlib;
  extern double PH;
  extern int CHARGES_PH_INDEPENDENT_FLAG, MAX_SEQ_POSITION;
  int i, j, k, q, n, m, d, nn, dd, numrot, first, flag;
  extern double **CHARGE_PRODUCT;

  first = 1;
  fixed_res_VARIABLE_POSITION = (VARIABLE_POSITION *)calloc(MAX_RESIDUES, sizeof(VARIABLE_POSITION));
  i=1; num_residues=0;
  while(protein->Template[i].seq_position!=ENDFLAG)
    {
      while(protein->Template[i].atom_ptr->atomnumber != 8 && protein->Template[i].seq_position!=ENDFLAG)
	++i;
      if(strcmp(protein->Template[i].residuetype, "PRO")!=0 && strcmp(protein->Template[i].residuetype, "GLY")!=0 && 
		protein->Template[i].seq_position!=ENDFLAG)
	{
	  ++num_residues;
	  fixed_res_VARIABLE_POSITION[num_residues].seq_position = protein->Template[i].seq_position;
	  while(protein->Template[i].seq_position == fixed_res_VARIABLE_POSITION[num_residues].seq_position && protein->Template[i].seq_position!=ENDFLAG)
	    ++i;
	}
      else
	while(protein->Template[i].seq_position == protein->Template[i-1].seq_position && protein->Template[i].seq_position!=ENDFLAG)
	  ++i;
    }
  fixed_res_VARIABLE_POSITION[num_residues+1].seq_position = ENDFLAG;
  
	/* filter out PRO/GLY/CYD */
	original_num_varpos=protein->parameters.numVarPositions;
  i=1; 
  while(protein->var_pos[i].seq_position!=ENDFLAG)
  {
	flag=0;
  	j=1;
  	while(fixed_res_VARIABLE_POSITION[j].seq_position!=ENDFLAG)
	{
		if(fixed_res_VARIABLE_POSITION[j].seq_position == protein->var_pos[i].seq_position)
			flag =1;
		++j;
	}
	if(flag == 0)
	{
		
		fprintf(stderr,"WARNING cannot find residue %s in Template structure or it is a PRO, GLY, or CYD....will ignore\n",
    protein->seqpos_text_map[seqpos_to_inputted_string(protein->var_pos[i].seq_position, protein->seqpos_text_map)].seqpos_text);			
		/*											
		exit(1);
		*/

		protein->var_pos[i].seq_position=MAX_SEQ_POSITION+protein->var_pos[i].seq_position;
		--protein->parameters.numVarPositions;
	}
	++i;
  }

   i = 1;
   sort_VARIABLE_POSITION(&first, &original_num_varpos, protein->var_pos);
   protein->var_pos[protein->parameters.numVarPositions+1].seq_position = ENDFLAG; /* mark the end of the array */


 /* extract native rotamers and calculate, internal energies, and link them to fixed_res_VARIABLE_POSITION */   

  pdbATOM_to_VARIABLE_POSITION(fixed_res_VARIABLE_POSITION, protein->Template, protein->resparam);
  fixed_res_RESPARAM = (RESPARAM *)calloc(MAX_RESIDUES, sizeof(RESPARAM));
  for(k=1;k<=num_residues;++k)
    fixed_res_RESPARAM[k] = *fixed_res_VARIABLE_POSITION[k].choice[1].resparam_ptr;
  strcpy(fixed_res_RESPARAM[num_residues + 1].residuetype, "zzz");
  intrinsic_rotamer_energy(fixed_res_RESPARAM);
  for(k=1;k<=num_residues;++k)
    *fixed_res_VARIABLE_POSITION[k].choice[1].resparam_ptr = fixed_res_RESPARAM[k];
  free_memory(fixed_res_RESPARAM);
    
  tempseqVarpos1 = (int *)calloc(MAX_RESIDUES, sizeof(int));
  tempseqVarpos = (int *)calloc(MAX_RESIDUES, sizeof(int));
  seqVarpos = (int *)calloc(MAX_RESIDUES, sizeof(int));
  

	/* user-defined variable positions */
  user_defined_varpos = (VARIABLE_POSITION *)calloc((protein->parameters.numVarPositions+2), sizeof(VARIABLE_POSITION));
  i=1;
  while(protein->var_pos[i].seq_position!=ENDFLAG)
    {
      user_defined_varpos[i] = protein->var_pos[i];
      user_defined_varpos[i].neighbor_level = ENDFLAG; /* user-defined position */
      seqVarpos[i] = user_defined_varpos[i].seq_position;
      ++i;
    }
  user_defined_varpos[i].seq_position=ENDFLAG; 
  seqVarpos[i]=ENDFLAG; 
  free_memory(protein->var_pos);
  
  /* define the neighboring residues near the mutation sites that will be allowed to float  */ 
  
  if(protein->parameters.neighbordist != PI_PLUS_LN10_PLUS_SQRT2_MINUS_E_1e5)
  {
  	protein->parameters.numVarPositions = find_neighbors(seqVarpos, protein->Template, tempseqVarpos1,&protein->parameters.neighbordist, 1.5); 
	k = protein->parameters.numVarPositions;

	/* second shell of neighbors */
  	protein->parameters.numVarPositions = find_neighbors(tempseqVarpos1, protein->Template, tempseqVarpos, &protein->parameters.neighbordist, 1.0);
   }
   else /* float only nearest-neighbors */
	{
		protein->parameters.neighbordist = 0;
		protein->parameters.numVarPositions = find_neighbors(seqVarpos, protein->Template, tempseqVarpos1, &protein->parameters.neighbordist, 1.5);

	k = protein->parameters.numVarPositions;

	for(i=1;i<=k;++i)
 		tempseqVarpos[i] = tempseqVarpos1[i];
	tempseqVarpos[k+1] = ENDFLAG;	


	}

  neighbor_varpos =  (VARIABLE_POSITION *)calloc((protein->parameters.numVarPositions+2), sizeof(VARIABLE_POSITION));
    
  listofvarpos_to_VARIABLE_POSITION(tempseqVarpos, neighbor_varpos, protein->Template, protein->parameters.numVarPositions, protein->resparam);

  protein->parameters.numVarPositions = num_residues; 
  protein->var_pos = (VARIABLE_POSITION *)calloc((protein->parameters.numVarPositions+2), sizeof(VARIABLE_POSITION));
    
	/* meld the three VARIABLE_POSITION arrays together into protein->var_pos */

  meld_VARIABLE_POSITION(protein->var_pos, fixed_res_VARIABLE_POSITION, neighbor_varpos, user_defined_varpos, fixed_positions);

  protein->numVarPositions =  protein->parameters.numVarPositions;


  /* mark level 1 neighbors */
  for(j=1;j<=k;++j)
  {
  	i=1;
  	while(protein->var_pos[i].seq_position!=tempseqVarpos1[j])
		++i;
	if(protein->var_pos[i].neighbor_level != ENDFLAG)
		protein->var_pos[i].neighbor_level = 1;

  }
   
  
  i=1;  q=1;
  while(protein->var_pos[i].seq_position!=ENDFLAG)
  {
	for(j=1;j<=protein->var_pos[i].number_of_choices;++j)
		protein->var_pos[i].choice[j].bkbn = &protein->var_pos[i].bkbn;

	while(protein->seqpos_text_map[q].seq_position!=protein->var_pos[i].seq_position)
		++q;
        protein->var_pos[i].seqpos_text_map_ptr = &protein->seqpos_text_map[q];	

	++i;
  }

  define_coresurfint(protein->Template, protein->var_pos);
   
  /* for ionizable residues, also include the neutral form as an option */
  if(CHARGES_PH_INDEPENDENT_FLAG == 1)
    { 
      i=1;  
      while(protein->var_pos[i].seq_position!=ENDFLAG)
	{
/*	  if(protein->var_pos[i].fixed_flag == 0)  */
	    {
	      fixed_res_RESPARAM = (RESPARAM *)calloc(10, sizeof(RESPARAM));
	      neutral_choice = (CHOICE *)calloc(10, sizeof(CHOICE)); 
	      k=0;
	      j=1;
	      while(j<=protein->var_pos[i].number_of_choices)
		{
		  if(protein->var_pos[i].choice[j].resparam_ptr->pKa!=20)
		    {
		      p=0;
		      if(PH != ENDFLAG)  /* if  PH = ENDFLAG, we'll deal w/ these energies elsewhere; just get both forms */
			{
			  if(protein->var_pos[i].choice[j].resparam_ptr->protonation_flag == 1)	/* protonated form was input */
			    p = exp(-(PH - protein->var_pos[i].choice[j].resparam_ptr->pKa)*LN10);
			  else	/* de-protonated form was input */
			    p = exp(-(protein->var_pos[i].choice[j].resparam_ptr->pKa - PH)*LN10);
			  p = p/(1.0 + p);	/* probab of this form at this pH */
			  p = 1.0 -p;  /* probab of other form at this pH */
			}
		      if(p>=0.001 || PH == ENDFLAG) /* If the other form is at least 0.1%, include it */
			{   
			  ++k;
			  neutral_choice[k].bkbn = &protein->var_pos[i].bkbn;
		      protein->var_pos[i].fixed_flag=0;
			  q=1;
			  while(protein->var_pos[i].choice[j].resparam_ptr->residuecode != -protein->resparam[q].residuecode)
			    ++q;
			  neutral_choice[k].resparam_ptr = (RESPARAM *)malloc(sizeof(RESPARAM));
			  *neutral_choice[k].resparam_ptr = protein->resparam[q];
		     
			  /*  neutral and charged variants should both have the same number of rotamers, unless the natural
			      rotamer was added on to the charged form in meld_VARIABLE_POSITION  */
			  if(protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->number_of_rotamers != 
			     neutral_choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers)
			    {
			      dummy_rotamerlib = *neutral_choice[k].resparam_ptr->rotamerlib_ptr;
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr = (ROTAMERLIB *)malloc(sizeof(ROTAMERLIB));
			      *neutral_choice[k].resparam_ptr->rotamerlib_ptr = dummy_rotamerlib;

			
				if(protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->number_of_rotamers != 1)
			    {    
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers = 
				neutral_choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers + 1;
			      numrot = neutral_choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers;
			    }
			    else
				{
					neutral_choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers = 1;
					numrot = neutral_choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers;
				}

			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->freq_array = (double *)calloc((numrot + 3), sizeof(double));
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer = (ROTAMER *)calloc((numrot + 3), sizeof(ROTAMER));
			    

				if(protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->number_of_rotamers != 1)
			    {
			      max_freq=0; denom=0;
			      for(n=1;n<=(numrot-1);++n)
				{
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n] = dummy_rotamerlib.rotamer[n];
				  if(dummy_rotamerlib.rotamer[n].freq > max_freq)
				    max_freq = dummy_rotamerlib.rotamer[n].freq;
				  denom += dummy_rotamerlib.rotamer[n].freq;	
				}
				}
				else
				{
					max_freq = 1.0;
					denom = 0;
				}
				

			
			      n=1;
			      while(protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->rotamer[n].native_rotamer_flag != 1)
				++n;

			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot] = 
				protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->rotamer[n];

			    
			      for(q=0;q<=neutral_choice[k].resparam_ptr->rotamerlib_ptr->numChi+1;++q)
				{
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].chi[q] = 
				    protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->rotamer[n].chi[q];	
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].std_dev[q] = 
				    protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->rotamer[n].std_dev[q];
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].half_std_dev[q] = 
				    0.49*protein->var_pos[i].choice[j].resparam_ptr->rotamerlib_ptr->rotamer[n].std_dev[q];
				}
			     
			      /* calculate rotamerlets */
			      if(neutral_choice[k].resparam_ptr->rotamerlib_ptr->num_rotamerlets != 0)
				{	
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet = (double **)calloc(neutral_choice[k].resparam_ptr->rotamerlib_ptr->num_rotamerlets+2, sizeof(double *));
				  for(q=1;q<=neutral_choice[k].resparam_ptr->rotamerlib_ptr->num_rotamerlets;++q)
				    {
				      neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q] = (double *)calloc(neutral_choice[k].resparam_ptr->rotamerlib_ptr->numChi+2, sizeof(double));

				      for(n=0;n<=neutral_choice[k].resparam_ptr->rotamerlib_ptr->numChi+1;++n)
					neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][n] = neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].chi[n];
				    }
	  
				  if(neutral_choice[k].resparam_ptr->rotamerlib_ptr->num_rotamerlets >= 3)
				    {
				      q=1;
				      for(m=-1;m<=1;++m)
					{
					  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][1] = 
					    neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].chi[1] + m*NATURAL_ROTAMER_SPREAD;
					  if(neutral_choice[k].resparam_ptr->rotamerlib_ptr->num_rotamerlets>=9)
					    {
					      n=q;
					      for(d=-1;d<=1;++d)
						{
						  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][1] = 
						    neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[n][1];
						  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][2] = 
						    neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][2] + d*NATURAL_ROTAMER_SPREAD;
			  
						  if(neutral_choice[k].resparam_ptr->rotamerlib_ptr->num_rotamerlets>=27)
						    {
						      nn = q;
						      for(dd = -1;dd<=1;++dd)
							{
							  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][1] = neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[nn][1];
							  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][2] = neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[nn][2];
							  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][3] = neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].rotamerlet[q][3] + dd*NATURAL_ROTAMER_SPREAD;
							  ++q;						    
							}
						    }
						  else	 					    
						    ++q;						    
						}						     
					    }
					  else
					    ++q;
					}
				    }
				}
			   
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[numrot].freq = max_freq;
			      denom += max_freq;
		
			      for(n=1;n<=numrot;++n)
				neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].freq = 
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].freq/denom;

			  

				if(neutral_choice[k].resparam_ptr->ligand_flag==0)
		     		sort_ROTAMER(&first,&numrot,neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer);
				else
				{
					if(neutral_choice[k].resparam_ptr->rotamerlib_ptr->numChi==1)
					{
						for(n=1;n<=numrot;++n)
							neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].chi[1] = n;
					}
				}




			
			      for(n=1;n<=numrot;++n)
				neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].index = n;
		
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->freq_array[0] = 0; 
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->freq_array[1] = 0;
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->freq_array[2] = 0;
		
			      for(n=1;n<=numrot;++n)
				{
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].rank = n;
				  neutral_choice[k].resparam_ptr->rotamerlib_ptr->freq_array[n] = 
				    neutral_choice[k].resparam_ptr->rotamerlib_ptr->rotamer[n].freq + 
				    neutral_choice[k].resparam_ptr->rotamerlib_ptr->freq_array[n-1];
				}
			      neutral_choice[k].resparam_ptr->rotamerlib_ptr->freq_array[numrot+1] = ENDFLAG;
			    }
			  fixed_res_RESPARAM[k] = *neutral_choice[k].resparam_ptr;	
			}
		    }
		  ++j;
		}
	      if(k!=0)	/*  at least one ionizable residue at this position, so re-adjust the composition and freq variables */
		{
		  strcpy(fixed_res_RESPARAM[k+1].residuetype, "zzz");	/*  get intrinsic rotamer energies  */
		  intrinsic_rotamer_energy(fixed_res_RESPARAM);
		  for(n=1;n<=k;++n)
		    {
		      *neutral_choice[n].resparam_ptr = fixed_res_RESPARAM[n];
		      protein->var_pos[i].choice[protein->var_pos[i].number_of_choices + n] = neutral_choice[n];
		    }    
		
		  protein->var_pos[i].number_of_choices += k;
		
		  /*  generate residue_freq array  */
		  protein->var_pos[i].residue_freqs[0]=0;
		  protein->var_pos[i].residue_freqs[1]=0;
		  protein->var_pos[i].residue_freqs[2]=0;
		  for(q=1;q<=protein->var_pos[i].number_of_choices;++q)
		    {   
		      protein->var_pos[i].choice[q].composition = 1.0/((double)protein->var_pos[i].number_of_choices);
		      protein->var_pos[i].residue_freqs[q] = protein->var_pos[i].residue_freqs[q-1] + protein->var_pos[i].choice[q].composition;
		    }
		  protein->var_pos[i].residue_freqs[protein->var_pos[i].number_of_choices+1] = ENDFLAG;
		}
		free_memory(neutral_choice);
	  	free_memory(fixed_res_RESPARAM);
	    }
	  ++i;
	  
	}
	
      if(PH!=ENDFLAG)
	{
	  /* termini partial atomic charges set to the charge state, pH dependent, but not location dependent */
	  k=1;
	  while(strcmp(protein->resparam[k].residuetype, "CTN")!=0)
	    ++k;
	  i=1;
	  while(strcmp(protein->resparam[i].residuetype,"CTE")!=0)
	    ++i;
	  p = exp(-(protein->resparam[i].pKa - PH)*LN10);
	  p = p/(1.0 + p);
	  for(q=1;q<MAX_RES_SIZE;++q)
	    {
	      protein->resparam[i].atom[q].charge = (1-p)*protein->resparam[k].atom[q].charge + p*protein->resparam[i].atom[q].charge;
	      protein->resparam[i].atom[q].charge2 = protein->resparam[i].atom[q].charge*protein->resparam[i].atom[q].charge;
	    }
	
	  k=1;
	  while(strcmp(protein->resparam[k].residuetype,"NTN")!=0)
	    ++k;
	  i=1;
	  while(strcmp(protein->resparam[i].residuetype,"NTE")!=0)
	    ++i;
	  p = exp(-(PH - protein->resparam[i].pKa)*LN10);
	  p = p/(1.0 + p);
	  for(q=1;q<MAX_RES_SIZE;++q)
	    {
	      protein->resparam[i].atom[q].charge = (1-p)*protein->resparam[k].atom[q].charge + p*protein->resparam[i].atom[q].charge;
	      protein->resparam[i].atom[q].charge2 = protein->resparam[i].atom[q].charge*protein->resparam[i].atom[q].charge;
	    }
	
	  i=1;
	  while(strcmp(protein->resparam[i].residuetype,"zzz")!=0)
	    {
	      for(n=4;n<=(protein->resparam[i].numberofAtoms+3);++n)
		{
		  j=1;
		  while(strcmp(protein->resparam[j].residuetype,"zzz")!=0)
		    {
		      for(k=4;k<=(protein->resparam[j].numberofAtoms+3);++k)
			{
			  CHARGE_PRODUCT[protein->resparam[i].atom[n].charge_class][protein->resparam[j].atom[k].charge_class] = 
			    protein->resparam[i].atom[n].charge*protein->resparam[j].atom[k].charge;	
			
			}   
		      ++j;
		    }
		}	
	      ++i;	
	    }
	}  /* end pH dependent terminii stuff */

    } /* end PH dependence stuff */
     
    
    
  /* set some pointers and defaults */
  i=1; q=1; protein->parameters.num_total_choices = 0;
  while(protein->var_pos[i].seq_position!=ENDFLAG)
    {
      protein->var_pos[i].varpos_index = i;
      protein->var_pos[i].lookup_energy_ptr = NULL;
      protein->var_pos[i].resimer = NULL;
      protein->var_pos[i].dee_index = ENDFLAG;
	
      
	protein->parameters.num_total_choices += protein->var_pos[i].number_of_choices;

	protein->var_pos[i].number_of_resimers = 0;  
      for(k=1;k<=protein->var_pos[i].number_of_choices;++k)
      {
	protein->var_pos[i].choice[k].lookup_res_ptr = NULL;
	protein->var_pos[i].number_of_resimers += protein->var_pos[i].choice[k].resparam_ptr->rotamerlib_ptr->number_of_rotamers;
	protein->var_pos[i].choice[k].in_use_flag = 1;
      }
      ++i;	
    }
    
  
	/* extract fixed atoms, calculate their energies */
	if(protein->fixed_atoms==NULL)
        {
	  protein->fixed_atoms = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM));
	  protein->mini_fixed_atoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
	  
	  /* Born radii from Template structure assigned to fixed atoms */
	  energy_calc(protein->mini_Template, 1, 0);

	  i=1;
	  while(protein->mini_Template[i].seq_position!=ENDFLAG)
	    {
	      protein->Template[i].born_radius = protein->mini_Template[i].born_radius;
	      protein->Template[i].sasa = 0;
	      ++i;
	    }
	  protein->Template[i].sasa = ENDFLAG;
	  protein->Template[i].seq_position = ENDFLAG;
	  protein->Template[i].atom_ptr= NULL;

	  make_fixed_pdbATOM(protein->Template, protein->var_pos, protein->fixed_atoms);
	
	  pdbATOM_to_mini_pdbATOM(protein->fixed_atoms, protein->mini_fixed_atoms);
	  
	  protein->fixedatoms_energy = fixedatoms_energy(protein->mini_fixed_atoms, protein->var_pos, protein->resparam);
	    
	  protein->parameters.fixedatoms_energy = ENERGY_TOTAL(protein->fixedatoms_energy);

		
	}

  free_memory(fixed_res_VARIABLE_POSITION); 
  free_memory(neighbor_varpos);
  free_memory(user_defined_varpos);
  free_memory(tempseqVarpos1);
  free_memory(tempseqVarpos);
  free_memory(seqVarpos);
}

/* similar to Marshall & Mayo JMB 305: 619-631 (2001) */
  /* Place a linear 3-carbon pseudosidechain (methyl acetylene MEA) at each position; get born radii */
  /* If a pseudosidechain has sasa <= CORESURF_MEA_SURF and born-CG >= CORESURF_MEA_BORN  call it core */
  /* surface if sasa < CORESURF_MEA_SURF and born-CG < CORESURF_MEA_BORN
	else interfacial */

void define_coresurfint(pdbATOM *Template, VARIABLE_POSITION *var_pos)
{
    mini_pdbATOM  *mea_minipdb;
    pdbATOM *mea_pdb, *sideatom;
    double *dummychi;
    int i,j,k,q,flag;
    double sasa;
    extern RESPARAM *MEA;
    extern double CORESURF_MEA_BORN, CORESURF_MEA_SASA;
    SIDECHAIN side;

    mea_minipdb = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
    mea_pdb = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM)); 
    sideatom = (pdbATOM *)calloc(MAX_RES_SIZE, sizeof(pdbATOM));
    dummychi = (double *)calloc(MEA->rotamerlib_ptr->numChi+2, sizeof(double));
    side.atom=(micro_pdbATOM *)calloc(MAX_RES_SIZE,sizeof(micro_pdbATOM));

    for(i=1;i<=MEA->rotamerlib_ptr->numChi;++i)
      dummychi[i] = 180;
    dummychi[0]=0; dummychi[MEA->rotamerlib_ptr->numChi+1]=0;     
    
    i=1; j=1; k=1;
    while(var_pos[i].seq_position!=ENDFLAG)
      {	
	flag=0;
	if(var_pos[i].choice[1].resparam_ptr->ligand_flag==0)
		flag=1;	
	if(var_pos[i].seqpos_text_map_ptr == NULL)
		flag=1;
	else
		if(strstr(var_pos[i].seqpos_text_map_ptr->seqpos_text,"l")==0)
			flag=1;
	
	if(flag==1)
	{
	while(Template[j].seq_position != var_pos[i].seq_position)
	  {
	    mea_pdb[k] = Template[j];
	    ++k; ++j;
	  }
	while(Template[j].seq_position == var_pos[i].seq_position)
	  {
	    if(Template[j].atom_ptr->other_info<0)
	      {
		mea_pdb[k] = Template[j];
		++k;
	      }	
	    ++j;
	  }

			
    	build_a_SideChain(&side, &var_pos[i].bkbn, MEA, &var_pos[i].seq_position, dummychi);
	make_side_pdbATOM(sideatom, &side);	
	q=1;
	while(sideatom[q].atom_ptr!=NULL)
	  {
	    mea_pdb[k] = sideatom[q];
	    ++q; ++k;
	  }

	}
	++i;
      }
    mea_pdb[k].seq_position=ENDFLAG; mea_pdb[k].sasa=ENDFLAG; mea_pdb[k].atom_ptr=NULL; strcpy(mea_pdb[k].residuetype, "END");
    pdbATOM_to_mini_pdbATOM(mea_pdb, mea_minipdb); 
    energy_calc(mea_minipdb, 1, 1);
   
   
    i=1; j=1; 
    while(var_pos[i].seq_position!=ENDFLAG)
      {
	
	flag=0;
	if(var_pos[i].choice[1].resparam_ptr->ligand_flag==0)
		flag=1;	
	if(var_pos[i].seqpos_text_map_ptr == NULL)
		flag=1;
	else
		if(strstr(var_pos[i].seqpos_text_map_ptr->seqpos_text,"l")==0)
			flag=1;
	
	if(flag==1)
	{    
	while(mea_minipdb[j].seq_position!=var_pos[i].seq_position)
	  ++j;
	k=j;    
	sasa = 0;
	while(mea_minipdb[j].seq_position == var_pos[i].seq_position)
	  {
	    if(mea_minipdb[j].atom_ptr->other_info>0 || mea_minipdb[j].atom_ptr->atomnumber ==8)
	      sasa += mea_minipdb[j].sasa;
	    ++j;
	  }
	
	    while(strcmp(mea_minipdb[k].atom_ptr->atomname, "CG")!=0)  
	      ++k; 
	    var_pos[i].mea_sasa = sasa; 
	    var_pos[i].mea_born =  mea_minipdb[k].born_radius;
	
	if(mea_minipdb[k].born_radius >= CORESURF_MEA_BORN && var_pos[i].mea_sasa <= CORESURF_MEA_SASA) //  core
		var_pos[i].core_flag = 'c';
	else if(mea_minipdb[k].born_radius < CORESURF_MEA_BORN && var_pos[i].mea_sasa > CORESURF_MEA_SASA) //  surface
		var_pos[i].core_flag = 's';
	else
		var_pos[i].core_flag = 'i';  	// interfacial

	  }
	else
		var_pos[i].core_flag = 'l';
	++i;
      }

	free_memory(mea_minipdb); free_memory(mea_pdb); free_memory(dummychi); 
	free_memory(sideatom); free_memory(side.atom);

}