/* EGAD: SASA.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 solvent-accessible-surface-area (SASA) calculations.
   
Many of these functions are adapted/taken from Scott LeGrand's sasa program 
	Le Grand, S. and Merz, K. "A Rapid Approximation to Molecular Surface Area via
        the Use of Boolean Logic and Look-Up Tables" J Comp Chem 14: 349 (1993)	
*/


#include "SASA.h"

SASA_SUM SASA_SUM_TEMP;

/* Matrix of contact distances squared by atom class */
double **SASA_RADII_SQRD_MATRIX=NULL;
double **SASA_RADII2_DIFF=NULL;	    /* SASA_RADII2_DIFF[i][j] = sasa_radius2_i - sasa_radius2_j */
CARTESIAN_VECTOR ***VECTOR_PTR=NULL;

double SASA_DISTANCE_CUTOFF;

int atm_pts_minus_MAXPTS; 
int  *BOX_OFFSET;		/* Non-bond box offsets */

/* SASA mask look-up tables */

unsigned char SASA_CLOSEST_POINT[100][200];

SASA_point *point; /* sphere points for sasa calculations */

double	ATOM_POINTS;		/* # of points on sphere */
int  sasa_count[258];		/* Binary->Decimal conv. table */

SASA_ATOM   **box;    /* Non-bond box */


/* reads the points file containing the SASA lookup table and initializes arrays */

void SASA_readfiles(char *SASA_points_file, ATOMPARAM *atomparam) 
{
FILE 		*fp;		/* Generic file pointer */
int 		fpnum;		/* Integer file pointer */
int	 	i, j, k,  l;		/* for loop counters */
int nj,nk; /* Used to set up non-bond box offsets for SASA calculations */
int		atom_resolution;	/* Grid resolution */
int atmpts, calc_points_flag;
char *filename;

if(SASA_RADII_SQRD_MATRIX!=NULL)
	return;

filename=NULL;

    box = (SASA_ATOM **)calloc(MAX_ATOMS, sizeof(SASA_ATOM *));

    BOX_OFFSET = (int *)calloc(27,sizeof(int));


    point = (SASA_point *)calloc(SASA_MAXPOINTS,sizeof(SASA_point));


    SASA_RADII_SQRD_MATRIX = (double **)calloc(MAXATOMTYPES, sizeof(double *));
    SASA_RADII2_DIFF = (double **)calloc(MAXATOMTYPES, sizeof(double *));
    for(i=0;i<MAXATOMTYPES;++i)
    {
	SASA_RADII_SQRD_MATRIX[i] = (double *)calloc(MAXATOMTYPES, sizeof(double));
	SASA_RADII2_DIFF[i] = (double *)calloc(MAXATOMTYPES, sizeof(double));
    }


atom_resolution=200;

SASA_DISTANCE_CUTOFF=0.0;

for(i=1;i<MAXATOMTYPES;++i)
{
    SASA_DISTANCE_CUTOFF=max(atomparam[i].sasa_radius,SASA_DISTANCE_CUTOFF);
    for(j=1;j<MAXATOMTYPES;++j)
	{
	    SASA_RADII2_DIFF[i][j] = atomparam[i].sasa_radius2 - atomparam[j].sasa_radius2;
	    SASA_RADII_SQRD_MATRIX[i][j]= pow((atomparam[i].sasa_radius + atomparam[j].sasa_radius), 2);
	}
}

/* Set sasa distance cutoff to 2xMaximum Atomic Radius */
  SASA_DISTANCE_CUTOFF = 2*SASA_DISTANCE_CUTOFF;

/* Determine non-bond box offsets */
  l=0;
  for (i= -1;i<=0;i++)
    {
      if (i==-1)
        nj=1;
      else
        nj=0;
      for(j= -1;j<=nj;j++)
        {
          if(i== -1 || j==-1)
            nk=1;
          else
            nk=0;
          for (k= -1;k<=nk;k++)
            {
              BOX_OFFSET[l]=i*1000000+j*1000+k;
              if (l!=13)
                BOX_OFFSET[l+14]= -BOX_OFFSET[l];
              l++;
            }
        }
    }   
    calc_points_flag = 0;

	if(SASA_points_file==NULL)
		calc_points_flag = 1;		
	else if(does_this_file_exist(SASA_points_file)==0)
		calc_points_flag = 1;	

	if(calc_points_flag == 1)
	{
		generate_SASA_points(SASA_points_file); 

	}
	else
	{
		fp=fopen_file(SASA_points_file,"r");
	
		fpnum=fileno(fp);
		read(fpnum,&atmpts,4);
		atm_pts_minus_MAXPTS = atmpts-SASA_MAXPOINTS;

  		read(fpnum,&atom_resolution,4);

  		for (i=0;i<atmpts;i++)
  			read(fpnum,point[i].mask,32*atom_resolution);
  		read(fpnum,SASA_CLOSEST_POINT,20000);
		read(fpnum,sasa_count,1024);
		fclose(fp);	

		ATOM_POINTS = (double)atmpts;
  	}

	
}




int nrnd(int num)		/* Generate random int from 0 to num-1 */
{
    double f ;
  f=(double)(num)*(double)(random())/(2147483648.0);
  return ((int)f) ;
}


int cmp_box(SASA_ATOM *a1, SASA_ATOM *a2)	/* Compares a1 and a2, returns TRUE if a1>a2 */
{
	if (a1->box>a2->box)
		return(TRUE);
	else
		if (a1->box<a2->box)
			return(FALSE);
		else
			if (a1->number>a2->number)
				return(TRUE);
			else
				return(FALSE);
}

void sort_box(int lower,int upper)	/* Quicksort routine */
{
int		pivot;		/* Quicksort pivet point */
SASA_ATOM	*pivot_atom;	/* Pointer to pivot atom */
SASA_ATOM	*temp;		/* Temporary atom pointer */
int		uc,lc;		/* Upper and lower counters */
int		flag;		/* Loop flag */

	/* Randomly Determine Pivot atom */
	  pivot=nrnd(upper-lower+1)+lower;  
	
	pivot_atom=box[pivot];
	/* Place pivot atom at top of list so it doesn't get swapped */
	temp=box[pivot];
	box[pivot]=box[upper];
	box[upper]=temp;
	/* Initialize counters */
	lc=lower;
	uc=upper;
	flag=TRUE;
	/* Perform swaps around pivot point */
	do
	{
		while (cmp_box(pivot_atom,box[lc]))
			lc++;
		while (cmp_box(box[uc],pivot_atom))
			uc--;
		if (lc<uc)
		{
			temp=box[uc];
			box[uc]=box[lc];
			box[lc]=temp;
		}
		else
			flag=FALSE;
	}
	while (flag);
	/* Now, recursively sort respective halves of box */
	if (uc-1>lower)
		sort_box(lower,uc-1);
	if (upper>uc)
		sort_box(uc,upper);
}


int find_box(int box_id,int *boxstart,int *boxend, int box_atoms)		/* Locate a non-bond box and its boundaries */
{
    int		upper;		/* Upper boundary of search */
    int		lower;		/* lower boundary of search */
    int		middle;		/* Midpoint of search */
    int		successful;	/* Flag indicating success of search */
    
	/* Initialize search */
	lower=0;
	upper=box_atoms-1;
	middle=(upper+lower)/2;
	successful=TRUE;
	/* Now conduct binary search for a box_id box */
	while (box[middle]->box!=box_id && successful)
	{
		if (upper<=lower)
			successful=FALSE;
		if (box[middle]->box<box_id)
			lower=middle+1;
		else
			upper=middle-1;
		middle=(upper+lower)/2;
	}
	if (successful)		/* Found desired box, now determine boundaries */
	{
		*boxstart=middle;
		*boxend=middle;
		while (*boxstart>=0 && box[*boxstart]->box==box_id)
			(*boxstart)--;
		(*boxstart)++;
		while (*boxend<box_atoms && box[*boxend]->box==box_id)
			(*boxend)++;
		return(TRUE);
	}
	else
		return(FALSE);
}


/* scales components in SASA_SUM by scale factor c */

SASA_SUM sasa_sum_scale(double c, SASA_SUM b)	 
{   
    SASA_SUM a;
    
    a.sasa_total = (b.sasa_total)*c;    
    a.sp2 = (b.sp2)*c;	 
    a.sp3_S = (b.sp3_S)*c;	 
    a.O = (b.O)*c;  
    a.N = (b.N)*c; 
    a.H = (b.H)*c; 
    
    return(a); 
}

/* calculates SASA_SUM for mini_pdbATOM array pdb; also calculates SASA for each atom --> pdb[i].sasa */
SASA_SUM SASAcalc(mini_pdbATOM *pdb)
{
int		i,j,k,l,m;	/* Counter variables */
int		boxstart;	/* Beginning of active non-bond box */
int		boxend;		/* End of active non-bond box */
int		box_id;		/* ID of active non-bond box */
int		cboxstart;	/* Beginning of temporary non-bond box */
int		cboxend;	/* End of temporary non-bond box */
int		cbox_id;	/* ID of temporary non-bond box */
SASA_ATOM	*a1,*a2;	/* Temporary atom pointers */
unsigned int	*p;		/* Pointer to look-up table masks */
double		phi,theta,r, r2;	/* Angular and distance variables */
double dx, dy, dz;
mini_pdbATOM *pdb0;
int		box_atoms;		/* Total atoms in area calculation */
double sasa_pts;
static SASA_ATOM *atm=NULL;		/* Atom array */
static int max_atoms=0;
SASA_SUM sasa_sum;
extern int BHY_CLASS;



	SASA_SUM_0(sasa_sum);
	
	if(atm == NULL)
	{
	    atm = (SASA_ATOM *)calloc(MAX_ATOMS, sizeof(SASA_ATOM));
	    max_atoms=MAX_ATOMS;
	}
	if(max_atoms<MAX_ATOMS) /* MAX_ATOMS has changed; resize static array */
	{
		free_memory(atm);

		atm = (SASA_ATOM *)calloc(MAX_ATOMS, sizeof(SASA_ATOM));
	    	max_atoms=MAX_ATOMS;
	}

	pdb0 = pdb;
	++pdb;
	box_atoms=0;
	while(pdb->seq_position!=ENDFLAG)
	{
	    pdb->sasa=0;
	    /* Only consider non-zero radius atoms */
	    if(pdb->atom_ptr->atom_ptr->sasa_radius>0.0)
	    {
		atm[box_atoms].pdb = pdb;
			
			/* Calculate non bond box coordinates */
			atm[box_atoms].box=
			((int)(atm[box_atoms].pdb->coord.x/SASA_DISTANCE_CUTOFF)+500)*1000000	
			+((int)(atm[box_atoms].pdb->coord.y/SASA_DISTANCE_CUTOFF)+500)*1000
			+(int)(atm[box_atoms].pdb->coord.z/SASA_DISTANCE_CUTOFF)+500
			-(atm[box_atoms].pdb->coord.x<0)*1000000
			-(atm[box_atoms].pdb->coord.y<0)*1000
			-(atm[box_atoms].pdb->coord.z<0);
			
			
			
		atm[box_atoms].number=box_atoms;
	    
		/* Initialize Atomic masks */
		atm[box_atoms].mask1=atm[box_atoms].mask2=atm[box_atoms].mask3=atm[box_atoms].mask4=
		atm[box_atoms].mask5=atm[box_atoms].mask6=atm[box_atoms].mask7=atm[box_atoms].mask8=0xffffffff;

			/* make non-bond box */
			box[box_atoms]= &atm[box_atoms];
			
		    ++box_atoms;	
	    }
	    ++pdb;
	}
	pdb = pdb0;

	
	if(box_atoms<2)
	{
	    if(box_atoms==1)  /* isolated atom; else no atom */
	    {
		atm[0].pdb->sasa = atm[0].pdb->atom_ptr->atom_ptr->sasa_sphere;
		sasa_sum.sasa_total = atm[0].pdb->sasa;
	    }		
	    return(sasa_sum);	
	}

	/* Now sort non-bond box */
	    sort_box(0,box_atoms-1);

	/* Do actual surface area calculation */
	for (boxstart=0;boxstart<box_atoms;)
	{
		/* Determine which box we're in */
		box_id=box[boxstart]->box;
		/* Now find boundary of active non-bond box */
		boxend=boxstart;
		while (boxend<box_atoms && box[boxend]->box==box_id)
			boxend++;
		/* Now calculate interactions with active box's */
		/* 13 neighbors and itself */
		for (l=0;l<14;l++)
		{
			cbox_id=box_id+BOX_OFFSET[l];
			/* Test if neighboring box is present */
			if (find_box(cbox_id,&cboxstart,&cboxend, box_atoms))
			{	/* Calculate cbox<->box interactions */
				for (i=boxstart;i<boxend;i++)
				{
					a1=box[i];
					/* If in same box, don't calculate self-self interactions */
					if (l==13)
						cboxend=i;
					for (j=cboxstart;j<cboxend;j++)
					{
						a2=box[j];
						/* Calculate distance between a1 and a2 */
						dx = a2->pdb->coord.x - a1->pdb->coord.x;
						dy = a2->pdb->coord.y - a1->pdb->coord.y;
						dz = a2->pdb->coord.z - a1->pdb->coord.z;
						r2 = dx*dx + dy*dy + dz*dz;
						/* Do these two atoms overlap? */
						if(r2>0)
						if(r2<SASA_RADII_SQRD_MATRIX[a1->pdb->atom_ptr->atom_ptr->atom_class][a2->pdb->atom_ptr->atom_ptr->atom_class])
						{
							/* determine distance and orientation */
							r=sqrt(r2);


							if(r + a1->pdb->atom_ptr->atom_ptr->sasa_radius >= a2->pdb->atom_ptr->atom_ptr->sasa_radius)
							if(r + a2->pdb->atom_ptr->atom_ptr->sasa_radius >= a1->pdb->atom_ptr->atom_ptr->sasa_radius)
							{
							dx/=r; dy/=r; dz/=r;
							theta=acos(dz);
							phi=atan2(dx,dy);
							if (phi<0.0)
								phi+=PI_2;

					/* a1 buried by a2 */
					/* we want the area of BHY atoms, yet do not want them to influence the area of other atoms */
							if(a2->pdb->atom_ptr->atom_ptr->atom_class != BHY_CLASS)
							{
							/* Now grab data from look-up table */
							k=SASA_CLOSEST_POINT[(int)(PI_under100*theta)][(int)(PI_under100*phi)];
							m=(int)(141.42136*sqrt(1.0-(SASA_RADII2_DIFF[a1->pdb->atom_ptr->atom_ptr->atom_class][a2->pdb->atom_ptr->atom_ptr->atom_class]+r2)/(a1->pdb->atom_ptr->atom_ptr->sasa_diameter*r)));
							/* Now do ANDing with a1's masks */
							p= &(point[k].mask[m][0]);
							a1->mask1=a1->mask1&(*(p++));
							a1->mask2=a1->mask2&(*(p++));
							a1->mask3=a1->mask3&(*(p++));
							a1->mask4=a1->mask4&(*(p++));
							a1->mask5=a1->mask5&(*(p++));
							a1->mask6=a1->mask6&(*(p++));
							a1->mask7=a1->mask7&(*(p++));
							a1->mask8=a1->mask8&(*(p++));
							}

							/* Now switch orientation to a2 */
							phi+=PI;
							if (phi>=PI_2)
								phi-=PI_2;
							theta=PI-theta;

				/* a2 buried by a1 */
			/* we want the area of BHY atoms, yet do not want them to influence the area of other atoms */
						if(a1->pdb->atom_ptr->atom_ptr->atom_class != BHY_CLASS)
							{
							/* Grab look-up table data for a2 */
						 	k=SASA_CLOSEST_POINT[(int)(PI_under100*theta)][(int)(PI_under100*phi)]; 
              						m=(int)(141.42136*sqrt(1.0-(SASA_RADII2_DIFF[a2->pdb->atom_ptr->atom_ptr->atom_class][a1->pdb->atom_ptr->atom_ptr->atom_class]+r2)/(a2->pdb->atom_ptr->atom_ptr->sasa_diameter*r)));
  							/* Do the same ANDing with a2's masks */
 							p= &(point[k].mask[m][0]);
 							a2->mask1=a2->mask1&(*(p++));
 							a2->mask2=a2->mask2&(*(p++));
 							a2->mask3=a2->mask3&(*(p++));
 							a2->mask4=a2->mask4&(*(p++));
 							a2->mask5=a2->mask5&(*(p++));
 							a2->mask6=a2->mask6&(*(p++));
 							a2->mask7=a2->mask7&(*(p++));
 							a2->mask8=a2->mask8&(*(p++));
							}

							}
						}
					}
				}
			}
		}
		boxstart=boxend;	
	}
	/* Now sum up invididual atomic surface areas */
	for (i=0;i<box_atoms;i++)
	{
			/* Count total number of exposed surface points */
			sasa_pts=sasa_count[atm[i].mask1&0x000000ff]+sasa_count[atm[i].mask1>>8&0x000000ff]+
			sasa_count[atm[i].mask1>>16&0x000000ff]+sasa_count[atm[i].mask1>>24&0x000000ff]+
			sasa_count[atm[i].mask2&0x000000ff]+sasa_count[atm[i].mask2>>8&0x000000ff]+
			sasa_count[atm[i].mask2>>16&0x000000ff]+sasa_count[atm[i].mask2>>24&0x000000ff]+
			sasa_count[atm[i].mask3&0x000000ff]+sasa_count[atm[i].mask3>>8&0x000000ff]+
			sasa_count[atm[i].mask3>>16&0x000000ff]+sasa_count[atm[i].mask3>>24&0x000000ff]+
			sasa_count[atm[i].mask4&0x000000ff]+sasa_count[atm[i].mask4>>8&0x000000ff]+
			sasa_count[atm[i].mask4>>16&0x000000ff]+sasa_count[atm[i].mask4>>24&0x000000ff]+
			sasa_count[atm[i].mask5&0x000000ff]+sasa_count[atm[i].mask5>>8&0x000000ff]+
			sasa_count[atm[i].mask5>>16&0x000000ff]+sasa_count[atm[i].mask5>>24&0x000000ff]+
			sasa_count[atm[i].mask6&0x000000ff]+sasa_count[atm[i].mask6>>8&0x000000ff]+
			sasa_count[atm[i].mask6>>16&0x000000ff]+sasa_count[atm[i].mask6>>24&0x000000ff]+
			sasa_count[atm[i].mask7&0x000000ff]+sasa_count[atm[i].mask7>>8&0x000000ff]+
			sasa_count[atm[i].mask7>>16&0x000000ff]+sasa_count[atm[i].mask7>>24&0x000000ff]+
			sasa_count[atm[i].mask8&0x000000ff]+sasa_count[atm[i].mask8>>8&0x000000ff]+
			sasa_count[atm[i].mask8>>16&0x000000ff]+sasa_count[atm[i].mask8>>24&0x000000ff]+
			(atm_pts_minus_MAXPTS);			
      			atm[i].pdb->sasa=atm[i].pdb->atom_ptr->atom_ptr->sasa_sphere*sasa_pts/ATOM_POINTS;
			
			sasa_sum.sasa_total += atm[i].pdb->sasa;

			if(atm[i].pdb->atom_ptr->atom_ptr->torsiontype[0] != 'U')
			{

			/*
				if(atm[i].atom_ptr->donorflag[0] == 'B' || atm[i].atom_ptr->donorflag[0] == 'A' ||
					atm[i].atom_ptr->donorflag[0] == 'D')
						if(atm[i].pdb->sasa>0)
							atm[i].pdb->hbond_satisfied = 1;
			*/


				sasa_sum.E_transfer += atm[i].pdb->sasa*atm[i].pdb->atom_ptr->atom_ptr->asp;

				
				if(atm[i].pdb->atom_ptr->atom_ptr->torsiontype[0]=='C')
				{
					if(atm[i].pdb->atom_ptr->atom_ptr->torsiontype[1]=='T')
						sasa_sum.sp3_S += atm[i].pdb->sasa;
					else
						sasa_sum.sp2 += atm[i].pdb->sasa;
				}
				else
					if(atm[i].pdb->atom_ptr->atom_ptr->torsiontype[0]=='S')
						sasa_sum.sp3_S += atm[i].pdb->sasa;
					else
						if(atm[i].pdb->atom_ptr->atom_ptr->torsiontype[0]=='N')
							sasa_sum.N += atm[i].pdb->sasa;
						else
							if(atm[i].pdb->atom_ptr->atom_ptr->torsiontype[0]=='O')
								sasa_sum.O += atm[i].pdb->sasa;
							else
								if(atm[i].pdb->atom_ptr->atom_ptr->torsiontype[0]=='H')
									sasa_sum.H += atm[i].pdb->sasa;

			}
	}		

	


	return(sasa_sum);
}


/* calculates CARTESIAN_VECTOR for atompairs from a pair of mini_pdbATOM arrays 
	use for pairwise area calcs (presently not implemented) */
void pairwise_minipdbATOM_to_CARTESIAN_VECTOR(mini_pdbATOM *pdb, mini_pdbATOM *pdb2, CARTESIAN_VECTOR *vector)
{
  CARTESIAN_VECTOR *dummy_vector;
  mini_pdbATOM *dummypdb, *dummypdb2;
  int heavy_ctr1, heavy_ctr2;
  double dx, dy, dz, r2, r;
  
  dummy_vector = vector;
  vector->next = NULL;
  
  dummypdb = pdb;
  dummypdb2 = pdb2;
  ++pdb; 
  
  heavy_ctr1=0; 
  while(pdb->seq_position!=ENDFLAG)
    {
      if(pdb->atom_ptr->atom_ptr->sasa_radius>0)
	{
	  ++heavy_ctr1;
	  pdb2 =  dummypdb2;
	  ++pdb2;	
	  heavy_ctr2=0;
	  while(pdb2->seq_position!=ENDFLAG)
	    {
	      if(pdb2->atom_ptr->atom_ptr->sasa_radius>0)
		{
		  ++heavy_ctr2;
		  dx = pdb2->coord.x - pdb->coord.x;
		  dy = pdb2->coord.y - pdb->coord.y;
		  dz = pdb2->coord.z - pdb->coord.z;
		  r2 = dx*dx + dy*dy + dz*dz;
	 	  if(r2<=SASA_RADII_SQRD_MATRIX[pdb->atom_ptr->atom_ptr->atom_class][pdb2->atom_ptr->atom_ptr->atom_class]) 
		   if(r2>0) 
		    {
		      
		      r = sqrt(r2);
			
				/* triangle inequality; pointed out by Mark Voorhies */
			  if(r + pdb->atom_ptr->atom_ptr->sasa_radius >= pdb2->atom_ptr->atom_ptr->sasa_radius)
				if(r + pdb2->atom_ptr->atom_ptr->sasa_radius >= pdb->atom_ptr->atom_ptr->sasa_radius)
							{	


			  vector->i = heavy_ctr1;
		      vector->j = heavy_ctr2;
		      vector->r2 = r2;
		      vector->r = r;
		      dx/=r; dy/=r; dz/=r;
		      vector->i_theta = acos(dz);
		      vector->i_phi = atan2(dx,dy);
		      
		      dx = -dx; dy= -dy; dz= -dz;
		      vector->j_theta = acos(dz);
		      vector->j_phi = atan2(dx,dy);
		      
		      vector->next = (CARTESIAN_VECTOR *)malloc(sizeof(CARTESIAN_VECTOR));
		      vector = vector->next;
		      vector->next = NULL;

				}
		    }		    	    	    
		}
	      ++pdb2;
	    }   
	}
      ++pdb;
    }
  
  vector = dummy_vector;  
  
  pdb = dummypdb;
  pdb2 =  dummypdb2;
}

/* calculates CARTESIAN_VECTOR for atompairs from a mini_pdbATOM array 
	use for the self part of pairwise area calcs (presently not implemented) */
void minipdbATOM_to_CARTESIAN_VECTOR(mini_pdbATOM *pdb, CARTESIAN_VECTOR *vector)
{
  CARTESIAN_VECTOR *dummy_vector;
  mini_pdbATOM *dummypdb, *pdb2;
  int heavy_ctr1, heavy_ctr2;
  double dx, dy, dz, r2, r;
  
  dummy_vector = vector;
  vector->next = NULL;
  
  dummypdb = pdb;
  ++pdb;
  
  heavy_ctr1=0;
  while(pdb->seq_position!=ENDFLAG)
    {
      if(pdb->atom_ptr->atom_ptr->sasa_radius>0)
	{
	  ++heavy_ctr1;
	  pdb->box = ((int)(pdb->coord.x/SASA_DISTANCE_CUTOFF)+500)*1000000	
	    +((int)(pdb->coord.y/SASA_DISTANCE_CUTOFF)+500)*1000
	    +(int)(pdb->coord.z/SASA_DISTANCE_CUTOFF)+500
	    -(pdb->coord.x<0)*1000000
	    -(pdb->coord.y<0)*1000
	    -(pdb->coord.z<0);
	  pdb2 = pdb+1; heavy_ctr2 = heavy_ctr1+1; 	
	  while(pdb2->seq_position!=ENDFLAG)
	    {
	      if(pdb2->atom_ptr->atom_ptr->sasa_radius>0)
		{
		  dx = pdb2->coord.x - pdb->coord.x;
		  dy = pdb2->coord.y - pdb->coord.y;
		  dz = pdb2->coord.z - pdb->coord.z;
		  r2 = dx*dx + dy*dy + dz*dz;
		  if(r2<=SASA_RADII_SQRD_MATRIX[pdb->atom_ptr->atom_ptr->atom_class][pdb2->atom_ptr->atom_ptr->atom_class]) 
		    {
		      r = sqrt(r2);
	
			/* triangle inequality; pointed out by Mark Voorhies */
			  if(r + pdb->atom_ptr->atom_ptr->sasa_radius >= pdb2->atom_ptr->atom_ptr->sasa_radius)
				if(r + pdb2->atom_ptr->atom_ptr->sasa_radius >= pdb->atom_ptr->atom_ptr->sasa_radius)
							{	

				
		      vector->i = heavy_ctr1;
		      vector->j = heavy_ctr2;
		      vector->r2 = r2;

		      vector->r = r;
		      dx/=r; dy/=r; dz/=r;
		      vector->i_theta = acos(dz);
		      vector->i_phi = atan2(dx,dy);
		      
		      dx = -dx; dy= -dy; dz= -dz;
		      vector->j_theta = acos(dz);
		      vector->j_phi = atan2(dx,dy);
		      
		      vector->next = (CARTESIAN_VECTOR *)malloc(sizeof(CARTESIAN_VECTOR));
		      vector = vector->next;
		      vector->next = NULL;

				}
		    }	
		  ++heavy_ctr2;		    	    	    
		}
	      ++pdb2;
	    }   
	}
      ++pdb;
    }
  
  vector = dummy_vector;  
  
  pdb = dummypdb;
}


 
/* calculates SASA_SUM for mini_pdbATOM array pdb_ptr[]; also calculates SASA for each atom.  
  uses a precalculated CARTESIAN_VECTOR structure (must be calculated a priori for pdb_ptr[] )   */
SASA_SUM precalcSASAcalc(mini_pdbATOM *pdb_ptr[])
{
  int		i,j,k,l,m;	/* Counter variables */
  int		boxstart;	/* Beginning of active non-bond box */
  int		boxend;		/* End of active non-bond box */
  int		box_id;		/* ID of active non-bond box */
  int		cboxstart;	/* Beginning of temporary non-bond box */
  int		cboxend;	/* End of temporary non-bond box */
  int		cbox_id;	/* ID of temporary non-bond box */
  SASA_ATOM	*a1,*a2;	/* Temporary atom pointers */
  unsigned int	*p;		/* Pointer to look-up table masks */
  double phi,theta;	/* Angular and distance variables */
  CARTESIAN_VECTOR *a1a2_vector;
  mini_pdbATOM *dummypdb;
  int		box_atoms;		/* Total atoms in area calculation */
  int  n;
  double sasa_pts;
  static SASA_ATOM *atm = NULL;		/* Atom array */
  static int max_atoms;
  SASA_SUM  sasa_sum;
  
	
  SASA_SUM_0(sasa_sum);
  
  if(atm == NULL)
  {
    atm = (SASA_ATOM *)calloc(MAX_ATOMS, sizeof(SASA_ATOM));
    max_atoms = MAX_ATOMS;
  }
  if(max_atoms<MAX_ATOMS)	/* MAX_ATOMS has gotten bigger; resize static arrays */
  {
	free_memory(atm);

	atm = (SASA_ATOM *)calloc(MAX_ATOMS, sizeof(SASA_ATOM));
    	max_atoms = MAX_ATOMS;
  }  



  box_atoms=0; 
 
     n=0; 
     while(pdb_ptr[n]!=NULL)
     {
	dummypdb = pdb_ptr[n];
      while(dummypdb->seq_position!=ENDFLAG)	/* get the first mini_pdbATOM array */
	{  
	    if(dummypdb->atom_ptr->atom_ptr->sasa_radius>0)
	    { 
	      atm[box_atoms].atom_ptr = dummypdb->atom_ptr->atom_ptr;
	      atm[box_atoms].box = dummypdb->box;
	      atm[box_atoms].pdb = dummypdb;
	      
	      /* Initialize Atomic masks */
	      atm[box_atoms].mask1=atm[box_atoms].mask2=atm[box_atoms].mask3=atm[box_atoms].mask4=
		atm[box_atoms].mask5=atm[box_atoms].mask6=atm[box_atoms].mask7=atm[box_atoms].mask8=0xffffffff;
	      
	      /* make non-bond box */
	      box[box_atoms]= &atm[box_atoms];
	      
	      atm[box_atoms].number =  box_atoms+1;  
	      ++box_atoms;
	     }     
	      ++dummypdb; 
	}
	++n;
    }	
  


  if(box_atoms<2)
    {
      if(box_atoms==1)  /* isolated atom; else no atom */
      {
	sasa_sum.sasa_total=atm[0].atom_ptr->sasa_sphere;
	atm[0].pdb->sasa = sasa_sum.sasa_total;
      }	
      return(sasa_sum);	
    }
  
  /* Now sort non-bond box */
  sort_box(0,box_atoms-1);
  


  /* Do actual surface area calculation */
  for (boxstart=0;boxstart<box_atoms;)
    {
      /* Determine which box we're in */
      box_id=box[boxstart]->box;
      /* Now find boundary of active non-bond box */
      boxend=boxstart;
      while (boxend<box_atoms && box[boxend]->box==box_id)
	boxend++;
      /* Now calculate interactions with active box's */
      /* 13 neighbors and itself */
      for (l=0;l<14;l++)
	{
	  cbox_id=box_id+BOX_OFFSET[l];
	  /* Test if neighboring box is present */
	  if(find_box(cbox_id,&cboxstart,&cboxend, box_atoms))
	    {	/* Calculate cbox<->box interactions */
	      for(i=boxstart;i<boxend;i++)
		{
		  a1=box[i];
		  /* If in same box, don't calculate self-self interactions */
		  if (l==13)
		    cboxend=i;
		  for (j=cboxstart;j<cboxend;j++)
		    {
	
		      a2=box[j];   

		      /* determine distance and orientation */ 
		      a1a2_vector = VECTOR_PTR[a1->number][a2->number];
		      if(a1a2_vector!=NULL)  /* NULL if they don't bury each other */
			{
					

			  if(a1->number<a2->number)  /* a1 is before a2; a1 is i in CARTESIAN_VECTOR */
			  {
			    theta=a1a2_vector->i_theta;
			    phi=a1a2_vector->i_phi;
			  }
			  else		/* a1 is after a2; a2 is j in CARTESIAN_VECTOR */
			  {
			    theta=a1a2_vector->j_theta;
			    phi=a1a2_vector->j_phi;
			  }  
			  if (phi<0.0)
			    phi+=PI_2;


			  /* Now grab data from look-up table */
			  k=SASA_CLOSEST_POINT[(int)(PI_under100*theta)][(int)(PI_under100*phi)];


			  m=(int)(141.42136*sqrt(1.0-(SASA_RADII2_DIFF[a1->atom_ptr->atom_class][a2->atom_ptr->atom_class]+a1a2_vector->r2)/(a1->atom_ptr->sasa_diameter*a1a2_vector->r)));
			  /* Now do ANDing with a1's masks */

			

			  p= &(point[k].mask[m][0]);
			  a1->mask1=a1->mask1&(*(p++));
			  a1->mask2=a1->mask2&(*(p++));
			  a1->mask3=a1->mask3&(*(p++));
			  a1->mask4=a1->mask4&(*(p++));
			  a1->mask5=a1->mask5&(*(p++));
			  a1->mask6=a1->mask6&(*(p++));
			  a1->mask7=a1->mask7&(*(p++));
			  a1->mask8=a1->mask8&(*(p++));


		

			  /* Now switch orientation to a2 */
			  phi+=PI;
			  if (phi>=PI_2)
			    phi-=PI_2;
			  theta=PI-theta;
			  /* Grab look-up table data for a2 */

		

			  k=SASA_CLOSEST_POINT[(int)(PI_under100*theta)][(int)(PI_under100*phi)]; 
			  m=(int)(141.42136*sqrt(1.0-(SASA_RADII2_DIFF[a2->atom_ptr->atom_class][a1->atom_ptr->atom_class]+a1a2_vector->r2)/(a2->atom_ptr->sasa_diameter*a1a2_vector->r)));
			  /* Do the same ANDing with a2's masks */
			  p= &(point[k].mask[m][0]);
			  a2->mask1=a2->mask1&(*(p++));
			  a2->mask2=a2->mask2&(*(p++));
			  a2->mask3=a2->mask3&(*(p++));
			  a2->mask4=a2->mask4&(*(p++));
			  a2->mask5=a2->mask5&(*(p++));
			  a2->mask6=a2->mask6&(*(p++));
			  a2->mask7=a2->mask7&(*(p++));
			  a2->mask8=a2->mask8&(*(p++)); 

		


			}
		    }
		}
	    }
	}
      boxstart=boxend;	
    }
  /* Now sum up invididual atomic surface areas */
  for (i=0;i<box_atoms;i++)
    {
      /* Count total number of exposed surface points */
      sasa_pts=sasa_count[atm[i].mask1&0x000000ff]+sasa_count[atm[i].mask1>>8&0x000000ff]+
	sasa_count[atm[i].mask1>>16&0x000000ff]+sasa_count[atm[i].mask1>>24&0x000000ff]+
	sasa_count[atm[i].mask2&0x000000ff]+sasa_count[atm[i].mask2>>8&0x000000ff]+
	sasa_count[atm[i].mask2>>16&0x000000ff]+sasa_count[atm[i].mask2>>24&0x000000ff]+
	sasa_count[atm[i].mask3&0x000000ff]+sasa_count[atm[i].mask3>>8&0x000000ff]+
	sasa_count[atm[i].mask3>>16&0x000000ff]+sasa_count[atm[i].mask3>>24&0x000000ff]+
	sasa_count[atm[i].mask4&0x000000ff]+sasa_count[atm[i].mask4>>8&0x000000ff]+
	sasa_count[atm[i].mask4>>16&0x000000ff]+sasa_count[atm[i].mask4>>24&0x000000ff]+
	sasa_count[atm[i].mask5&0x000000ff]+sasa_count[atm[i].mask5>>8&0x000000ff]+
	sasa_count[atm[i].mask5>>16&0x000000ff]+sasa_count[atm[i].mask5>>24&0x000000ff]+
	sasa_count[atm[i].mask6&0x000000ff]+sasa_count[atm[i].mask6>>8&0x000000ff]+
	sasa_count[atm[i].mask6>>16&0x000000ff]+sasa_count[atm[i].mask6>>24&0x000000ff]+
	sasa_count[atm[i].mask7&0x000000ff]+sasa_count[atm[i].mask7>>8&0x000000ff]+
	sasa_count[atm[i].mask7>>16&0x000000ff]+sasa_count[atm[i].mask7>>24&0x000000ff]+
	sasa_count[atm[i].mask8&0x000000ff]+sasa_count[atm[i].mask8>>8&0x000000ff]+
	sasa_count[atm[i].mask8>>16&0x000000ff]+sasa_count[atm[i].mask8>>24&0x000000ff]+
	(atm_pts_minus_MAXPTS);			
        atm[i].pdb->sasa = atm[i].atom_ptr->sasa_sphere*sasa_pts/ATOM_POINTS;
	sasa_sum.E_transfer += atm[i].pdb->sasa*atm[i].atom_ptr->asp;
	sasa_sum.sasa_total += atm[i].pdb->sasa;

	
	
			if(atm[i].atom_ptr->torsiontype[0] != 'U')
			{

		/*
				if(atm[i].atom_ptr->donorflag[0] == 'B' || atm[i].atom_ptr->donorflag[0] == 'A' ||
					atm[i].atom_ptr->donorflag[0] == 'D')
						if(atm[i].pdb->sasa>0)
							atm[i].pdb->hbond_satisfied = 1;
			*/

				if(atm[i].atom_ptr->torsiontype[0]=='C')
				{
					if(atm[i].atom_ptr->torsiontype[1]=='T')
						sasa_sum.sp3_S += atm[i].pdb->sasa;
					else
						sasa_sum.sp2 += atm[i].pdb->sasa;
				}
				else
					if(atm[i].atom_ptr->torsiontype[0]=='S')
						sasa_sum.sp3_S += atm[i].pdb->sasa;
					else
						if(atm[i].atom_ptr->torsiontype[0]=='N')
							sasa_sum.N += atm[i].pdb->sasa;
						else
							if(atm[i].atom_ptr->torsiontype[0]=='O')
								sasa_sum.O += atm[i].pdb->sasa;
							else
								if(atm[i].atom_ptr->torsiontype[0]=='H')
									sasa_sum.H += atm[i].pdb->sasa;

			}
    }		
    
	

    return(sasa_sum);
}

/* for a mini_pdbATOM array that already has sasa's calculated for each atom, return SASA_SUM */

SASA_SUM get_SASA_SUM(mini_pdbATOM *pdb)
{
	SASA_SUM sasa_sum;
	int i;

	SASA_SUM_0(sasa_sum);

	i=1;
	while(pdb[i].seq_position!=ENDFLAG)
	{
		sasa_sum.E_transfer += pdb[i].sasa*pdb[i].atom_ptr->atom_ptr->asp;
		sasa_sum.sasa_total += pdb[i].sasa;

			if(pdb[i].atom_ptr->atom_ptr->torsiontype[0] != 'U')
			{
				if(pdb[i].atom_ptr->atom_ptr->torsiontype[0]=='C')
				{
					if(pdb[i].atom_ptr->atom_ptr->torsiontype[1]=='T')
						sasa_sum.sp3_S += pdb[i].sasa;
					else
						sasa_sum.sp2 += pdb[i].sasa;
				}
				else
					if(pdb[i].atom_ptr->atom_ptr->torsiontype[0]=='S')
						sasa_sum.sp3_S += pdb[i].sasa;
					else
						if(pdb[i].atom_ptr->atom_ptr->torsiontype[0]=='N')
							sasa_sum.N += pdb[i].sasa;
						else
							if(pdb[i].atom_ptr->atom_ptr->torsiontype[0]=='O')
								sasa_sum.O += pdb[i].sasa;
							else
								if(pdb[i].atom_ptr->atom_ptr->torsiontype[0]=='H')
									sasa_sum.H += pdb[i].sasa;
			}

		++i;
	}

	return(sasa_sum);
}


/* free some memory */
void freeCARTESIAN_VECTOR(CARTESIAN_VECTOR *vector)
{
    CARTESIAN_VECTOR *temp;
    
    while(vector!=NULL)  
    { 
	temp=vector->next; 
	free_memory(vector);     
	vector=temp;      
    }
}

/* free some memory */
void free_CARTESIAN_VECTOR_ptr(CARTESIAN_VECTOR ***vector_ptr, int sizeof_array)
{
    int i;
    
    for(i=0;i<=sizeof_array;++i)
	{
	    if(vector_ptr[i]!=NULL)
		free_memory(vector_ptr[i]);
	    vector_ptr[i]=NULL;
	}
	
    free_memory(vector_ptr);
    vector_ptr=NULL;
    
}




/* SASA mask generation structure */
typedef struct {
  int distance,id;
  } SASA_DISTANCE_TYPE;

SASA_DISTANCE_TYPE *dmat[SASA_MAXPOINTS];

int cmp_dmat(SASA_DISTANCE_TYPE *d1, SASA_DISTANCE_TYPE *d2)	/* Compares two dmat entries, returns TRUE if d1>d2 */
{
	/* Perform Comparison */
	if (d1->distance>d2->distance)
		return(TRUE);
	else
		if (d1->distance<d2->distance)
			return(FALSE);
		else
			if (d1->id>d2->id)
				return(TRUE);
			else
				return(FALSE);
}


void sort_dmat(int lower, int upper)		/* Quicksort dmat entries by distance */
{
int		pivot;		/* Pivot point */
SASA_DISTANCE_TYPE	*pivot_entry;	/* Pointer to pivot entry */
SASA_DISTANCE_TYPE	*temp;		/* Temporary pointer for swaps */
int		uc,lc;		/* Upper and lower counters */
int		flag;		/* Loop flag */
	/* Randomly determine pivot point */
	pivot=nrnd(upper-lower+1)+lower;
	pivot_entry=dmat[pivot];
	/* Place pivot enry at top of list so it can't get swapped */
	temp=dmat[pivot];
	dmat[pivot]=dmat[upper];
	dmat[upper]=temp;
	/* Initialize Counters */
	lc=lower;
	uc=upper;
	flag=TRUE;
	/* Perform Swaps around pivot point */
	do
	{
		while (cmp_dmat(pivot_entry,dmat[lc]))
			lc++;
		while (cmp_dmat(dmat[uc],pivot_entry))
			uc--;
		if (lc<uc)
		{
			temp=dmat[uc];
			dmat[uc]=dmat[lc];
			dmat[lc]=temp;
		}
		else
			flag=FALSE;
	}
	while (flag);
	/* Now, recursively sort respective halves of dmat list */
	if (uc-1>lower)
		sort_dmat(lower,uc-1);
	if (upper>uc)
		sort_dmat(uc,upper);
}



void generate_SASA_points(char *SASA_points_file)
{
int		i,j,k,l;	/* Counter variables */
unsigned int mask[8],dmask;  /* temporary SASA masks */
double 		x,y,z;		/* Coordinates of current point */
double		x1,y1,z1;	/* Temporary point coordinates */
double		radius;		/* distance from current point */
double		dummy_radius;	/* Temporary distance variable */
double		phi,theta;	/* Phi and theta for closest point table */
int		idist;		/* scaled and rounded distance from current point */
int		old_idist;	/* old idist */
int		closest;	/* Temporary closest point variable */
FILE		*fp;		/* File pointer */
int		fpnum,atom_points;		/* Integer file pointer */
int atom_resolution = 200;

	/* define sphere points; 256 equidistance points on a sphere surface */
	atom_points = 255;

	point[0].x = -0.997867; point[0].y = 0.034699; point[0].z = -0.055298;
	point[1].x = -0.174403; point[1].y = 0.158603; point[1].z = 0.971817;
	point[2].x = -0.069099; point[2].y = -0.917882; point[2].z = 0.390792;
	point[3].x = -0.771924; point[3].y = -0.325910; point[3].z = -0.545817;
	point[4].x = -0.470000; point[4].y = 0.353700; point[4].z = 0.808700;
	point[5].x = 0.509116; point[5].y = 0.634720; point[5].z = -0.581318;
	point[6].x = -0.934269; point[6].y = -0.192894; point[6].z = 0.299890;
	point[7].x = 0.424412; point[7].y = -0.685120; point[7].z = -0.592017;
	point[8].x = 0.444213; point[8].y = -0.303009; point[8].z = 0.843125;
	point[9].x = -0.964695; point[9].y = -0.009500; point[9].z = -0.263199;
	point[10].x = -0.604336; point[10].y = 0.203712; point[10].z = -0.770246;
	point[11].x = 0.283395; point[11].y = 0.736686; point[11].z = 0.613988;
	point[12].x = -0.988599; point[12].y = -0.113100; point[12].z = 0.099400;
	point[13].x = 0.162603; point[13].y = -0.615210; point[13].z = 0.771412;
	point[14].x = -0.574391; point[14].y = -0.466293; point[14].z = -0.672790;
	point[15].x = 0.620482; point[15].y = 0.008200; point[15].z = 0.784178;
	point[16].x = -0.189890; point[16].y = 0.942451; point[16].z = -0.275186;
	point[17].x = -0.200089; point[17].y = 0.431277; point[17].z = -0.879753;
	point[18].x = -0.732259; point[18].y = -0.009399; point[18].z = 0.680962;
	point[19].x = -0.064801; point[19].y = 0.423409; point[19].z = 0.903618;
	point[20].x = 0.724400; point[20].y = 0.216200; point[20].z = -0.654600;
	point[21].x = -0.243494; point[21].y = 0.941178; point[21].z = 0.234295;
	point[22].x = -0.388913; point[22].y = -0.304410; point[22].z = 0.869529;
	point[23].x = 0.755037; point[23].y = 0.654732; point[23].z = 0.035302;
	point[24].x = -0.255911; point[24].y = 0.224210; point[24].z = -0.940340;
	point[25].x = 0.032600; point[25].y = 0.467497; point[25].z = -0.883394;
	point[26].x = -0.920022; point[26].y = 0.335608; point[26].z = 0.202305;
	point[27].x = 0.069099; point[27].y = 0.917882; point[27].z = -0.390792;
	point[28].x = -0.797814; point[28].y = -0.466108; point[28].z = -0.382407;
	point[29].x = 0.600390; point[29].y = -0.789687; point[29].z = 0.126198;
	point[30].x = -0.183398; point[30].y = -0.401095; point[30].z = -0.897490;
	point[31].x = 0.645408; point[31].y = -0.759710; point[31].z = -0.079301;
	point[32].x = 0.217694; point[32].y = 0.971175; point[32].z = -0.097098;
	point[33].x = 0.943726; point[33].y = -0.235007; point[33].z = 0.232707;
	point[34].x = -0.769246; point[34].y = 0.635238; point[34].z = -0.068804;
	point[35].x = -0.467419; point[35].y = -0.850435; point[35].z = 0.241410;
	point[36].x = -0.085199; point[36].y = -0.646793; point[36].z = 0.757892;
	point[37].x = 0.931238; point[37].y = -0.047402; point[37].z = -0.361315;
	point[38].x = -0.860095; point[38].y = 0.455797; point[38].z = -0.229099;
	point[39].x = -0.427284; point[39].y = 0.903766; point[39].z = 0.025199;
	point[40].x = 0.255911; point[40].y = -0.224210; point[40].z = 0.940340;
	point[41].x = -0.570727; point[41].y = -0.211110; point[41].z = 0.793538;
	point[42].x = -0.469916; point[42].y = 0.782127; point[42].z = 0.409214;
	point[43].x = -0.121195; point[43].y = -0.879462; point[43].z = -0.460280;
	point[44].x = 0.630071; point[44].y = -0.606072; point[44].z = -0.485477;
	point[45].x = 0.395906; point[45].y = 0.896713; point[45].z = 0.197903;
	point[46].x = -0.629694; point[46].y = -0.768493; point[46].z = 0.113599;
	point[47].x = -0.755037; point[47].y = -0.654732; point[47].z = -0.035302;
	point[48].x = -0.543112; point[48].y = 0.505011; point[48].z = 0.670815;
	point[49].x = 0.939866; point[49].y = 0.320988; point[49].z = -0.116696;
	point[50].x = -0.719375; point[50].y = -0.648077; point[50].z = -0.249991;
	point[51].x = 0.604336; point[51].y = -0.203712; point[51].z = 0.770246;
	point[52].x = 0.769197; point[52].y = -0.635297; point[52].z = 0.068800;
	point[53].x = 0.969885; point[53].y = 0.208197; point[53].z = 0.126398;
	point[54].x = 0.654810; point[54].y = 0.553908; point[54].z = 0.514208;
	point[55].x = 0.786886; point[55].y = -0.110698; point[55].z = 0.607089;
	point[56].x = 0.381212; point[56].y = -0.494016; point[56].z = 0.781425;
	point[57].x = 0.964695; point[57].y = 0.009500; point[57].z = 0.263199;
	point[58].x = -0.582205; point[58].y = 0.463904; point[58].z = -0.667706;
	point[59].x = -0.908000; point[59].y = -0.276700; point[59].z = -0.314600;
	point[60].x = -0.820984; point[60].y = 0.393292; point[60].z = -0.413892;
	point[61].x = 0.582205; point[61].y = -0.463904; point[61].z = 0.667706;
	point[62].x = 0.692194; point[62].y = 0.528395; point[62].z = -0.491596;
	point[63].x = 0.860095; point[63].y = -0.455797; point[63].z = 0.229099;
	point[64].x = 0.643328; point[64].y = -0.243511; point[64].z = -0.725832;
	point[65].x = -0.217694; point[65].y = -0.971175; point[65].z = 0.097098;
	point[66].x = 0.249311; point[66].y = -0.967142; point[66].z = 0.049802;
	point[67].x = 0.997867; point[67].y = -0.034699; point[67].z = 0.055298;
	point[68].x = 0.178896; point[68].y = 0.826284; point[68].z = -0.534089;
	point[69].x = 0.100203; point[69].y = -0.603918; point[69].z = -0.790723;
	point[70].x = -0.304907; point[70].y = 0.695415; point[70].z = -0.650714;
	point[71].x = -0.875320; point[71].y = 0.251306; point[71].z = 0.413110;
	point[72].x = -0.397596; point[72].y = -0.094399; point[72].z = 0.912692;
	point[73].x = 0.064801; point[73].y = -0.423409; point[73].z = -0.903618;
	point[74].x = -0.546823; point[74].y = 0.757832; point[74].z = -0.355915;
	point[75].x = -0.045399; point[75].y = 0.245296; point[75].z = -0.968385;
	point[76].x = 0.889029; point[76].y = -0.149905; point[76].z = 0.432614;
	point[77].x = 0.458404; point[77].y = 0.166502; point[77].z = 0.873008;
	point[78].x = 0.732259; point[78].y = 0.009399; point[78].z = -0.680962;
	point[79].x = -0.142701; point[79].y = -0.579604; point[79].z = -0.802306;
	point[80].x = 0.661193; point[80].y = -0.692193; point[80].z = -0.289297;
	point[81].x = -0.969885; point[81].y = -0.208197; point[81].z = -0.126398;
	point[82].x = 0.964487; point[82].y = -0.262097; point[82].z = 0.032700;
	point[83].x = -0.259891; point[83].y = 0.855171; point[83].z = 0.448485;
	point[84].x = -0.158096; point[84].y = 0.854280; point[84].z = -0.495188;
	point[85].x = 0.484094; point[85].y = -0.669592; point[85].z = 0.563293;
	point[86].x = 0.470000; point[86].y = -0.353700; point[86].z = -0.808700;
	point[87].x = -0.253084; point[87].y = -0.918744; point[87].z = 0.303081;
	point[88].x = -0.726983; point[88].y = 0.317993; point[88].z = -0.608586;
	point[89].x = -0.692194; point[89].y = -0.528395; point[89].z = 0.491596;
	point[90].x = 0.800780; point[90].y = -0.467688; point[90].z = -0.374191;
	point[91].x = -0.931238; point[91].y = 0.047402; point[91].z = 0.361315;
	point[92].x = -0.225186; point[92].y = 0.711657; point[92].z = 0.665459;
	point[93].x = 0.662607; point[93].y = -0.565306; point[93].z = 0.491305;
	point[94].x = -0.567708; point[94].y = 0.051601; point[94].z = 0.821611;
	point[95].x = 0.401295; point[95].y = 0.915390; point[95].z = -0.032000;
	point[96].x = 0.225392; point[96].y = 0.002100; point[96].z = 0.974266;
	point[97].x = 0.562505; point[97].y = 0.759207; point[97].z = 0.327403;
	point[98].x = 0.849321; point[98].y = 0.455311; point[98].z = -0.267107;
	point[99].x = 0.360123; point[99].y = 0.833554; point[99].z = 0.418927;
	point[100].x = 0.079098; point[100].y = -0.840577; point[100].z = -0.535886;
	point[101].x = 0.298000; point[101].y = 0.681399; point[101].z = -0.668499;
	point[102].x = 0.085199; point[102].y = 0.646793; point[102].z = -0.757892;
	point[103].x = -0.424387; point[103].y = 0.685080; point[103].z = 0.592082;
	point[104].x = 0.543112; point[104].y = -0.505011; point[104].z = -0.670815;
	point[105].x = 0.920022; point[105].y = -0.335608; point[105].z = -0.202305;
	point[106].x = 0.719221; point[106].y = -0.396812; point[106].z = -0.570317;
	point[107].x = -0.939866; point[107].y = -0.320988; point[107].z = 0.116696;
	point[108].x = 0.225202; point[108].y = -0.711607; point[108].z = -0.665507;
	point[109].x = -0.605703; point[109].y = -0.712903; point[109].z = 0.353402;
	point[110].x = 0.174403; point[110].y = -0.158603; point[110].z = -0.971817;
	point[111].x = -0.875413; point[111].y = -0.452307; point[111].z = -0.170502;
	point[112].x = -0.661193; point[112].y = 0.692193; point[112].z = 0.289297;
	point[113].x = 0.891028; point[113].y = 0.100803; point[113].z = 0.442614;
	point[114].x = -0.078499; point[114].y = 0.994883; point[114].z = 0.063599;
	point[115].x = 0.567708; point[115].y = -0.051601; point[115].z = -0.821611;
	point[116].x = 0.183398; point[116].y = 0.401095; point[116].z = 0.897490;
	point[117].x = 0.469994; point[117].y = -0.782090; point[117].z = -0.409195;
	point[118].x = 0.346394; point[118].y = -0.829086; point[118].z = 0.438893;
	point[119].x = 0.059500; point[119].y = 0.211800; point[119].z = 0.975500;
	point[120].x = -0.257100; point[120].y = -0.507500; point[120].z = 0.822401;
	point[121].x = -0.974896; point[121].y = 0.122099; point[121].z = 0.186199;
	point[122].x = -0.720362; point[122].y = 0.625467; point[122].z = -0.299784;
	point[123].x = -0.318300; point[123].y = 0.532501; point[123].z = 0.784301;
	point[124].x = -0.442294; point[124].y = 0.868189; point[124].z = 0.224997;
	point[125].x = -0.387187; point[125].y = 0.151495; point[125].z = 0.909470;
	point[126].x = 0.158096; point[126].y = -0.854280; point[126].z = 0.495188;
	point[127].x = 0.639925; point[127].y = 0.240209; point[127].z = 0.729928;
	point[128].x = 0.886672; point[128].y = -0.462385; point[128].z = -0.003500;
	point[129].x = -0.122796; point[129].y = -0.989170; point[129].z = -0.080398;
	point[130].x = -0.862927; point[130].y = -0.103503; point[130].z = 0.494615;
	point[131].x = -0.484094; point[131].y = 0.669592; point[131].z = -0.563293;
	point[132].x = -0.197803; point[132].y = -0.945515; point[132].z = -0.258604;
	point[133].x = 0.189890; point[133].y = -0.942451; point[133].z = 0.275186;
	point[134].x = 0.252991; point[134].y = 0.918767; point[134].z = -0.303089;
	point[135].x = -0.012100; point[135].y = 0.982286; point[135].z = -0.186997;
	point[136].x = 0.988599; point[136].y = 0.113100; point[136].z = -0.099400;
	point[137].x = -0.628419; point[137].y = -0.419413; point[137].z = 0.655120;
	point[138].x = 0.498690; point[138].y = 0.678186; point[138].z = 0.539789;
	point[139].x = -0.795914; point[139].y = 0.159903; point[139].z = 0.583911;
	point[140].x = 0.394308; point[140].y = 0.792117; point[140].z = -0.465910;
	point[141].x = -0.059500; point[141].y = -0.211800; point[141].z = -0.975500;
	point[142].x = 0.629694; point[142].y = 0.768493; point[142].z = -0.113599;
	point[143].x = 0.797814; point[143].y = 0.466108; point[143].z = 0.382407;
	point[144].x = -0.072501; point[144].y = -0.744205; point[144].z = -0.664005;
	point[145].x = -0.346394; point[145].y = 0.829086; point[145].z = -0.438893;
	point[146].x = -0.381212; point[146].y = 0.494016; point[146].z = -0.781425;
	point[147].x = 0.045399; point[147].y = -0.245296; point[147].z = 0.968385;
	point[148].x = -0.178896; point[148].y = -0.826284; point[148].z = 0.534089;
	point[149].x = -0.639925; point[149].y = -0.240209; point[149].z = -0.729928;
	point[150].x = -0.562505; point[150].y = -0.759207; point[150].z = -0.327403;
	point[151].x = -0.100203; point[151].y = 0.603918; point[151].z = 0.790723;
	point[152].x = 0.570727; point[152].y = 0.211110; point[152].z = -0.793538;
	point[153].x = -0.401295; point[153].y = -0.915390; point[153].z = 0.032000;
	point[154].x = 0.771899; point[154].y = 0.325999; point[154].z = 0.545799;
	point[155].x = 0.243494; point[155].y = -0.941178; point[155].z = -0.234295;
	point[156].x = 0.072501; point[156].y = 0.744205; point[156].z = 0.664005;
	point[157].x = 0.908000; point[157].y = 0.276700; point[157].z = 0.314600;
	point[158].x = 0.574391; point[158].y = 0.466293; point[158].z = 0.672790;
	point[159].x = 0.257100; point[159].y = 0.507500; point[159].z = -0.822401;
	point[160].x = 0.874615; point[160].y = 0.484208; point[160].z = -0.024300;
	point[161].x = 0.779219; point[161].y = 0.096902; point[161].z = 0.619215;
	point[162].x = 0.397596; point[162].y = 0.094399; point[162].z = -0.912692;
	point[163].x = -0.886672; point[163].y = 0.462385; point[163].z = 0.003500;
	point[164].x = -0.443702; point[164].y = -0.389502; point[164].z = -0.807104;
	point[165].x = -0.643328; point[165].y = 0.243511; point[165].z = 0.725832;
	point[166].x = -0.654810; point[166].y = -0.553908; point[166].z = -0.514208;
	point[167].x = -0.800780; point[167].y = 0.467688; point[167].z = 0.374191;
	point[168].x = -0.509087; point[168].y = -0.634683; point[168].z = 0.581385;
	point[169].x = 0.318300; point[169].y = -0.532501; point[169].z = -0.784301;
	point[170].x = -0.249405; point[170].y = 0.967118; point[170].z = -0.049801;
	point[171].x = 0.581432; point[171].y = 0.808545; point[171].z = 0.090505;
	point[172].x = -0.079098; point[172].y = 0.840577; point[172].z = 0.535886;
	point[173].x = -0.225392; point[173].y = -0.002100; point[173].z = -0.974266;
	point[174].x = -0.943726; point[174].y = 0.235007; point[174].z = -0.232707;
	point[175].x = 0.443702; point[175].y = 0.389502; point[175].z = 0.807104;
	point[176].x = -0.298000; point[176].y = -0.681399; point[176].z = 0.668499;
	point[177].x = 0.304907; point[177].y = -0.695415; point[177].z = 0.650714;
	point[178].x = 0.605703; point[178].y = 0.712903; point[178].z = -0.353402;
	point[179].x = -0.800988; point[179].y = 0.573591; point[179].z = 0.171497;
	point[180].x = 0.122796; point[180].y = 0.989170; point[180].z = 0.080398;
	point[181].x = 0.078499; point[181].y = -0.994883; point[181].z = -0.063599;
	point[182].x = -0.786886; point[182].y = 0.110698; point[182].z = -0.607089;
	point[183].x = 0.121096; point[183].y = 0.879473; point[183].z = 0.460286;
	point[184].x = -0.444213; point[184].y = 0.303009; point[184].z = -0.843125;
	point[185].x = 0.795868; point[185].y = -0.159894; point[185].z = -0.583976;
	point[186].x = -0.150304; point[186].y = -0.283808; point[186].z = 0.947028;
	point[187].x = 0.142701; point[187].y = 0.579604; point[187].z = 0.802306;
	point[188].x = 0.388913; point[188].y = 0.304410; point[188].z = -0.869529;
	point[189].x = -0.747382; point[189].y = -0.618285; point[189].z = 0.243194;
	point[190].x = -0.360108; point[190].y = -0.833519; point[190].z = -0.419009;
	point[191].x = 0.628419; point[191].y = 0.419413; point[191].z = -0.655120;
	point[192].x = 0.875320; point[192].y = -0.251306; point[192].z = -0.413110;
	point[193].x = 0.270612; point[193].y = -0.341615; point[193].z = -0.900038;
	point[194].x = 0.820984; point[194].y = -0.393292; point[194].z = 0.413892;
	point[195].x = -0.443209; point[195].y = -0.479309; point[195].z = 0.757515;
	point[196].x = 0.442294; point[196].y = -0.868189; point[196].z = -0.224997;
	point[197].x = 0.546823; point[197].y = -0.757832; point[197].z = 0.355915;
	point[198].x = -0.818191; point[198].y = -0.338696; point[198].z = 0.464595;
	point[199].x = -0.407899; point[199].y = 0.889299; point[199].z = -0.206800;
	point[200].x = -0.032600; point[200].y = -0.467497; point[200].z = 0.883394;
	point[201].x = -0.394393; point[201].y = -0.792086; point[201].z = 0.465891;
	point[202].x = -0.874615; point[202].y = -0.484208; point[202].z = 0.024300;
	point[203].x = -0.395906; point[203].y = -0.896713; point[203].z = -0.197903;
	point[204].x = -0.888991; point[204].y = 0.149898; point[204].z = -0.432696;
	point[205].x = -0.630071; point[205].y = 0.606072; point[205].z = 0.485477;
	point[206].x = -0.662607; point[206].y = 0.565306; point[206].z = -0.491305;
	point[207].x = 0.875413; point[207].y = 0.452307; point[207].z = 0.170502;
	point[208].x = 0.000000; point[208].y = 0.000000; point[208].z = 1.000000;
	point[209].x = 0.194099; point[209].y = 0.083800; point[209].z = -0.977396;
	point[210].x = 0.387102; point[210].y = -0.151501; point[210].z = -0.909505;
	point[211].x = -0.581432; point[211].y = -0.808545; point[211].z = -0.090505;
	point[212].x = 0.467380; point[212].y = 0.850463; point[212].z = -0.241390;
	point[213].x = 0.747382; point[213].y = 0.618285; point[213].z = -0.243194;
	point[214].x = -0.600390; point[214].y = 0.789687; point[214].z = -0.126198;
	point[215].x = -0.724400; point[215].y = -0.216200; point[215].z = 0.654600;
	point[216].x = 0.407899; point[216].y = -0.889299; point[216].z = 0.206800;
	point[217].x = 0.443209; point[217].y = 0.479309; point[217].z = -0.757515;
	point[218].x = -0.283395; point[218].y = -0.736686; point[218].z = -0.613988;
	point[219].x = 0.270801; point[219].y = 0.225701; point[219].z = 0.935802;
	point[220].x = -0.719221; point[220].y = 0.396812; point[220].z = 0.570317;
	point[221].x = 0.934269; point[221].y = 0.192894; point[221].z = -0.299890;
	point[222].x = 0.041799; point[222].y = -0.949982; point[222].z = -0.309494;
	point[223].x = 0.974896; point[223].y = -0.122099; point[223].z = -0.186199;
	point[224].x = -0.645408; point[224].y = 0.759710; point[224].z = 0.079301;
	point[225].x = -0.018500; point[225].y = 0.776788; point[225].z = -0.629490;
	point[226].x = 0.862927; point[226].y = 0.103503; point[226].z = -0.494615;
	point[227].x = 0.719327; point[227].y = 0.648124; point[227].z = 0.250009;
	point[228].x = 0.720407; point[228].y = -0.625406; point[228].z = 0.299803;
	point[229].x = 0.727006; point[229].y = -0.317903; point[229].z = 0.608605;
	point[230].x = -0.498690; point[230].y = -0.678186; point[230].z = -0.539789;
	point[231].x = 0.427284; point[231].y = -0.903766; point[231].z = -0.025199;
	point[232].x = -0.620482; point[232].y = -0.008200; point[232].z = -0.784178;
	point[233].x = 0.357695; point[233].y = 0.575292; point[233].z = 0.735590;
	point[234].x = 0.801034; point[234].y = -0.573524; point[234].z = -0.171507;
	point[235].x = 0.000000; point[235].y = 0.000000; point[235].z = -1.000000;
	point[236].x = 0.012100; point[236].y = -0.982286; point[236].z = 0.186997;
	point[237].x = 0.259903; point[237].y = -0.855210; point[237].z = -0.448405;
	point[238].x = -0.458404; point[238].y = -0.166502; point[238].z = -0.873008;
	point[239].x = -0.270807; point[239].y = -0.225606; point[239].z = -0.935824;
	point[240].x = 0.818191; point[240].y = 0.338696; point[240].z = -0.464595;
	point[241].x = -0.270612; point[241].y = 0.341615; point[241].z = 0.900038;
	point[242].x = -0.964487; point[242].y = 0.262097; point[242].z = -0.032700;
	point[243].x = -0.438919; point[243].y = 0.051702; point[243].z = -0.897038;
	point[244].x = 0.018501; point[244].y = -0.776837; point[244].z = 0.629430;
	point[245].x = -0.891049; point[245].y = -0.100794; point[245].z = -0.442575;
	point[246].x = -0.779219; point[246].y = -0.096902; point[246].z = -0.619215;
	point[247].x = -0.357695; point[247].y = -0.575292; point[247].z = -0.735590;
	point[248].x = -0.041799; point[248].y = 0.949982; point[248].z = 0.309494;
	point[249].x = 0.150304; point[249].y = 0.283808; point[249].z = -0.947028;
	point[250].x = 0.438919; point[250].y = -0.051702; point[250].z = 0.897038;
	point[251].x = -0.194099; point[251].y = -0.083800; point[251].z = 0.977396;
	point[252].x = -0.849321; point[252].y = -0.455311; point[252].z = 0.267107;
	point[253].x = -0.162603; point[253].y = 0.615210; point[253].z = -0.771412;
	point[254].x = 0.197803; point[254].y = 0.945515; point[254].z = 0.258604;
	point[255].x = 0.200107; point[255].y = -0.431315; point[255].z = 0.879730;

	for (i=0;i<atom_points;i++)
		dmat[i]=(SASA_DISTANCE_TYPE *)malloc(sizeof(SASA_DISTANCE_TYPE));

	for (i=0;i<atom_points;i++)	/* Loop over all sphere points */
	{
		for (j=0;j<8;j++) /* Initialize all points active */
   			mask[j]=0xffffffff;
		x=point[i].x;     /* Reference from current point */
		y=point[i].y;
		z=point[i].z;
		for (j=0;j<atom_points;j++)  /* Generate distances of each point from current */
		{
			x1=point[j].x;
			y1=point[j].y;
			z1=point[j].z;
			
			radius=sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1));
			

			dmat[j]->distance=(int)(radius*100.0+0.5);
			/* Makes sure distances are legal, fudge a bit if not */
			if (dmat[j]->distance==200)
				dmat[j]->distance=199;
			dmat[j]->id=j;
		}
    		sort_dmat(0,atom_points-1);  /* Sort matrix by distance */
		/* Now turn off bits in masks in the order of their distance */
		/* from the current point */
		/* Set distance covered to zero */
		old_idist=idist=0;
		for (j=0;j<atom_points;j++)	/* Cycle through sorted distance list */
		{
			do
			{	/* Turn off each point between old_idist and idist */
				/* from current point */
				k=dmat[j]->id;
				dmask=0x00000001<<(k%32);
				k/=32;
				mask[k]=mask[k]^dmask;
				j++;
			}
      			while (j<atom_points && dmat[j]->distance<=idist);
			if (j==atom_points)	/* Determine new upper distance limit */
			{
				idist=200;
        			for (k=0;k<8;k++)
					mask[k]=0;
      			}
     			else
			{
				idist=dmat[j]->distance;
				j--;
			}
			for (l=old_idist;l<idist;l++)  /* Change masks out to upper distance limit */
				for (k=0;k<8;k++)
					point[i].mask[l][k]=mask[k];
			old_idist=l;		/* Copy into old upper distance limit */
		}

	}
	/* Now determine closest point table */
	for (i=0;i<100;i++) 
	{
		for (j=0;j<200;j++)
		{
			/* Determine coordinates of each point in closest point matrix */
			phi=PI*(double)(j)/100.0;
			theta=PI*(double)(i)/100.0;
			x=cos(phi)*sin(theta);
			y=sin(phi)*sin(theta);
			z=cos(theta);
			radius=100.0;
			closest= -1;
			for (k=0;k<atom_points;k++)	/* Now find nearest neighbor */
			{
				x1=point[k].x;
				y1=point[k].y;
				z1=point[k].z;
				dummy_radius=sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1));
				if (dummy_radius<radius)
				{
					radius=dummy_radius;
					closest=k;
				}
			}
			SASA_CLOSEST_POINT[i][j]=closest;
		}
	}
	for (i=1;i<129;i*=2)  /* Map mask values to numerical equivalents */
		for (k=i;k<256;k+=i*2)
			for (j=k;j<k+i;j++)
				sasa_count[j]++;

	atm_pts_minus_MAXPTS = atom_points-SASA_MAXPOINTS;
	ATOM_POINTS = (double)atom_points;

	if(SASA_points_file!=NULL)
	{
		fp=fopen(SASA_points_file,"w");
		fpnum=fileno(fp);
		/* Write size parameters */
		write(fpnum,&atom_points,4);
 		write(fpnum,&atom_resolution,4);  

		/* Write masks */
		for (i=0;i<atom_points;i++)
			write(fpnum,point[i].mask,32*atom_resolution);
		/* Write closest point matrix */
 	 	write(fpnum,SASA_CLOSEST_POINT,20000);
		/* Write bit counting table */
  		write(fpnum,sasa_count,1024);
  		fclose(fp); 
	}

}