/* EGAD: energy_functions.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 calculating energies.

*/


#include "energy_functions.h"

double FORCEFIELD_DISTANCE_CUTOFF=10;	/* max distance for an interaction  */
double FORCEFIELD_DISTANCE_CUTOFF_SQRD=100;	/* FORCEFIELD_DISTANCE_CUTOFF^2 */

SUPER_ENERGY super_energy;


/* This function pre-calculates the components of desolvation due to atoms whose relative distances will 
	not change during any of the calculations in EGAD (atoms <=1,3 and atoms on the same rigid ring). 
	Called by read_forcefield.cpp: read_forcefield
 GB implementation of Dominy and Brooks J Phys Chem B 103:3765-3773 (1999) */

void born_bonded(RESPARAM *inputresparam)
{
  int j,i,q, n;
  double r, r2, r4;
  double a, b, *chi;
  BACKBONE bkbn;
  SIDECHAIN side;
  pdbATOM *side_pdb;
  mini_pdbATOM *all_pdb;
  RESPARAM *resparam;
  int flag, bondstate;
  extern double BORN_P1;
  extern double BORN_P2;
  extern double BORN_P3;
  extern double BORN_P4;
  extern double ONE_OVER_BORN_LAMBDA;
  extern BACKBONE DUMMY_BKBN, PROLINE_BKBN;
  extern double *STATIONARY_LIGAND_TRANSFORM_MATRIX;

  n=1;
  while(strcmp(inputresparam[n].residuetype, "zzz")!=0)
  {
  resparam = &inputresparam[n];
  
  for(j=4;j<=(resparam->numberofAtoms+3);++j)
	resparam->atom[j].BORN_P4_volume = BORN_P4*resparam->atom[j].volume_bonded;
  
  
  chi = (double *)calloc(resparam->rotamerlib_ptr->numChi+2, sizeof(double));
  side_pdb = (pdbATOM *)calloc(MAX_RES_SIZE, sizeof(pdbATOM));
  all_pdb = (mini_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(mini_pdbATOM));
  side.atom = (micro_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(micro_pdbATOM));
  
  flag=0;
  if(resparam->residuetype[1]=='T' && resparam->residuetype[2]=='E'   /* strstr(resparam->residuetype, "TE")!=0 */)   /* terminii */
    flag=1;
    
    
  bkbn = DUMMY_BKBN;

  if(strcmp(resparam->residuetype,"PRO")==0)
    bkbn = PROLINE_BKBN;
  
/* set some aribitrary conformation in order to generate coordinates. 
	These really don't matter, since we are considering only <=1,3 interaction */

  chi[0] = 0; chi[resparam->rotamerlib_ptr->numChi+1] = 0; 
  for(q=1;q<=resparam->rotamerlib_ptr->numChi;++q)
    chi[q] = 180.0;
  

if(resparam->ligand_flag==1)
  {
	if(resparam->rotamerlib_ptr->numChi==1)	/* ligand "rotamers" are coords */
	{
		chi[1]=1;	/* just use the first ligand coord in the "rotamer" lib */
	}	
	else	/* ligand "rotamers" are transform matricies - use the base structure */
	    for(q=1;q<=12;++q)	
    		chi[q] = STATIONARY_LIGAND_TRANSFORM_MATRIX[q];
  }

	/* generate coordinates */
  q=1; 
  if(flag==0)
  {
    build_a_SideChain(&side, &bkbn,resparam,&q,chi);
    make_side_pdbATOM(side_pdb,&side);  
  }
  else
    strcpy(side_pdb[1].residuetype,"END");
  

  q=4;
  all_pdb[q].coord = bkbn.N; all_pdb[q].atom_ptr = &resparam->atom[q];
  all_pdb[q].sasa=0; all_pdb[q].seq_position=1;
  ++q;
  all_pdb[q].sasa=0; all_pdb[q].seq_position=1;
  all_pdb[q].coord = bkbn.H; all_pdb[q].atom_ptr = &resparam->atom[q];
  ++q;
  all_pdb[q].sasa=0; all_pdb[q].seq_position=1;
  all_pdb[q].coord = bkbn.CA; all_pdb[q].atom_ptr = &resparam->atom[q];
  ++q;
  all_pdb[q].sasa=0; all_pdb[q].seq_position=1;
  all_pdb[q].coord = bkbn.HA; all_pdb[q].atom_ptr = &resparam->atom[q];
  ++q;
  all_pdb[q].sasa=0; all_pdb[q].seq_position=1;
  all_pdb[q].coord = bkbn.CB; all_pdb[q].atom_ptr = &resparam->atom[q];
  ++q;
  i=1;
  while(strcmp(side_pdb[i].residuetype,"END")!=0)
    {
      all_pdb[q].coord = side_pdb[i].coord;
      all_pdb[q].atom_ptr = side_pdb[i].atom_ptr;
      all_pdb[q].sasa = 0; all_pdb[q].seq_position=1;
      ++q; ++i;
    }
  all_pdb[q].sasa=0; all_pdb[q].seq_position=1;
  all_pdb[q].coord = bkbn.C; all_pdb[q].atom_ptr = &resparam->atom[q];
  ++q;
  all_pdb[q].sasa=0; all_pdb[q].seq_position=1;
  all_pdb[q].coord = bkbn.O; all_pdb[q].atom_ptr = &resparam->atom[q];
  ++q;
  all_pdb[q].sasa=ENDFLAG; all_pdb[q].seq_position=ENDFLAG; all_pdb[q].atom_ptr=NULL;

  for(i=4;i<=(resparam->numberofAtoms+3);++i)
    if(resparam->atom[i].charge!=0)
	if(resparam->atom[i].atom_ptr->atom_class!=ENDFLAG)   
		resparam->atom[i].born_bonded = HALF_COULOMB_CONST*(1/resparam->atom[i].atom_ptr->united_radius)*(BORN_P1/resparam->atom[i].atom_ptr->united_radius - ONE_OVER_BORN_LAMBDA);
	
	/* loop over all atom-pairs in a residue (including backbone) */	
  r=0; bondstate=0;
  for(i=4;i<=(resparam->numberofAtoms+3);++i)
    if(resparam->atom[i].atom_ptr->united_radius>0)   /* only non-zero radius atoms */
      for(j=(i+1);j<=(resparam->numberofAtoms+3);++j) 
	if(resparam->atom[j].atom_ptr->united_radius>0)  /* only non-zero radius atoms */
	  {				
		/* find the connectivity */
	
	    /* consider here only atoms that are within <=2 bonds or within the same aromatic group */ 
	    if(resparam->atom[i].other_info>0 && resparam->atom[j].other_info>0 &&
	       resparam->atom[i].other_info == resparam->atom[j].other_info)  /* for fixed rings in a sidechain */
	      {
		r = Distance(all_pdb[i].coord, all_pdb[j].coord);
		if(resparam->atom[i].contactAtomNumber==j || resparam->atom[j].contactAtomNumber==i ||
		   resparam->atom[i].ringContactAtom==j || resparam->atom[j].ringContactAtom==i)
		    bondstate=1;
		else
		    if(resparam->atom[i].ringTwoBondAtom==j || resparam->atom[j].ringTwoBondAtom==i ||
		       resparam->atom[i].contactAtomNumber==resparam->atom[j].contactAtomNumber ||
		       resparam->atom[i].ringContactAtom==resparam->atom[j].ringContactAtom ||
		       resparam->atom[i].twobondAtomNumber==j || resparam->atom[j].twobondAtomNumber==i)
			bondstate=2;
		    else
			bondstate=3;   
	      }
	    else
	      if(strcmp(resparam->residuetype, "NTE")==0 &&   /* special case N-term HN-N */
		 ((strcmp(resparam->atom[i].atomname, "N")==0 && strcmp(resparam->atom[j].atomtype, "HN")==0) ||
		  (strcmp(resparam->atom[j].atomname, "N")==0 && strcmp(resparam->atom[i].atomtype, "HN")==0)))
		  {
			r = 1.01; bondstate=1;
		  }  	
	      else
		if(strcmp(resparam->residuetype, "CTE")==0 &&   /* special case C-term C-O */
		   ((strcmp(resparam->atom[i].atomname, "C")==0 && strncmp(resparam->atom[j].atomname, "O", 1)==0) ||
		    (strcmp(resparam->atom[j].atomname, "C")==0 && strncmp(resparam->atom[i].atomname, "O", 1)==0)))  
		  {  
			r = 1.229; bondstate=1;
		    }	
		 else
			if(strcmp(resparam->residuetype, "NTE")==0 &&   /* special case for the N-terminal group HN-CA interaction  */
			   ((strcmp(resparam->atom[i].atomname, "CA")==0 && strcmp(resparam->atom[j].atomtype, "HN")==0) || 
			    (strcmp(resparam->atom[j].atomname, "CA")==0 && strcmp(resparam->atom[i].atomtype, "HN")==0)))
			    {
				r = 2.09;  bondstate=2;   
			    }	 
			else      /* special case for the N-terminal HN-HN interaction */
			  if(strcmp(resparam->residuetype, "NTE")==0 && 
			     ((strcmp(resparam->atom[i].atomtype, "HN")==0 && strcmp(resparam->atom[j].atomtype, "HN")==0) &&
			      (resparam->atom[i].atomnumber!=resparam->atom[j].atomnumber)))
			    { 
				r = 1.7;  bondstate=2;  
			    }	 
			  else
			    if(strcmp(resparam->residuetype, "CTE")==0 &&   /* C-term CA-O */
			       ((strcmp(resparam->atom[i].atomname, "CA")==0 || strncmp(resparam->atom[j].atomname, "O", 1)==0) ||
				(strncmp(resparam->atom[i].atomtype, "O", 1)==0 || strcmp(resparam->atom[j].atomname, "CA")==0)))
			     { 
				r = 2.36;  bondstate=2; 
			     }	 
			    else
			      if(strcmp(resparam->residuetype, "CTE")==0 &&    /* C-term O-O */
				 ((strncmp(resparam->atom[i].atomname, "O", 1)==0 && strncmp(resparam->atom[j].atomtype, "O", 1)==0) &&
				  (resparam->atom[i].atomnumber!=resparam->atom[j].atomnumber)))
				  {
					r = 2.15;  bondstate=2; 
				  }   
		else		   
		  if(resparam->atom[i].contactAtomNumber==j)	
		  {    
		    r=resparam->atom[i].bondlength; bondstate=1;
		  }  
		  else
		    if(resparam->atom[j].contactAtomNumber==i)
		    {
		      r=resparam->atom[j].bondlength;	bondstate=1; 
		    }     
		    else
		      if(resparam->atom[i].ringContactAtom==j || resparam->atom[j].ringContactAtom==i ||
			 resparam->atom[i].ringTwoBondAtom==j || resparam->atom[j].ringTwoBondAtom==i)
			 {
				r = Distance(all_pdb[i].coord, all_pdb[j].coord); 
				if(resparam->atom[i].ringContactAtom==j || resparam->atom[j].ringContactAtom==i)
				    bondstate=1; 
				else
				     bondstate=2;    
			 }  	
			      else
				if(resparam->atom[i].twobondAtomNumber==j) 
				  {
				    a = resparam->atom[i].bondlength;
				    b = resparam->atom[resparam->atom[i].contactAtomNumber].bondlength;
				    r = sqrt(a*a + b*b - 2*a*b*cosd(resparam->atom[i].bondangle)); 
				    bondstate=2; 
				  }
				else
				  if(resparam->atom[j].twobondAtomNumber==i)    
				    {
				      a = resparam->atom[j].bondlength;
				      b = resparam->atom[resparam->atom[j].contactAtomNumber].bondlength;
				      r = sqrt(a*a + b*b - 2*a*b*cosd(resparam->atom[j].bondangle)); 
				      bondstate=2; 
				    }
				  else
				    if(resparam->atom[i].contactAtomNumber==resparam->atom[j].contactAtomNumber)
				      {	
					bondstate=2; 
					if(resparam->atom[i].other_info>0 || resparam->atom[j].other_info>0)
					  r = Distance(all_pdb[i].coord, all_pdb[j].coord); 
					else    /* Ni - Oi-1 */	
					  r = 2.2;   
				      }
				      else
				      if(resparam->atom[i].other_info == resparam->atom[j].other_info)  /* for fixed rings */
					{
					  r = Distance(all_pdb[i].coord, all_pdb[j].coord); 
					  bondstate=3; 
					}

		/* calculate desolvation based on the connectivity */
	    if(r!=0 && bondstate!=0)		
	    {
		r2=r*r; r4=r2*r2;
		if(resparam->atom[i].charge!=0 && strcmp(resparam->atom[j].atomname, "HU")!=0)
		{
		    if(bondstate==3)
			resparam->atom[i].born_bonded += born_nonbond_calc(&resparam->atom[i], &resparam->atom[j], r2);
		    else
			if(bondstate==1)
			    resparam->atom[i].born_bonded += BORN_P2*resparam->atom[j].volume_bonded/(r4);
			else
			    if(bondstate==2)    
				resparam->atom[i].born_bonded += BORN_P3*resparam->atom[j].volume_bonded/(r4);
		}
		if(resparam->atom[j].charge!=0 && strcmp(resparam->atom[i].atomname, "HU")!=0)
		{
		    if(bondstate==3)
			resparam->atom[j].born_bonded += born_nonbond_calc(&resparam->atom[j], &resparam->atom[i], r2);
		    else
			if(bondstate==1)
			    resparam->atom[j].born_bonded += BORN_P2*resparam->atom[i].volume_bonded/(r4);
			else
			    if(bondstate==2)    
				resparam->atom[j].born_bonded += BORN_P3*resparam->atom[i].volume_bonded/(r4);
		}
	        r=0; bondstate=0;
	    }
	}

	free_memory(chi);
	free_memory(side_pdb);
	free_memory(all_pdb);
	free_memory(side.atom);
	++n;
    }	  
}


/*
   This function precalculates the intrinsic energy for all rotamers.
   vdw energies between atoms of the same aromatic ring system are NOT calculated.
   1,4 and bonded born factors are added and stored in rotamer->born_factor array.
   Electrostatics between charged atoms are measured and stored in COULOMBIC list 
   called by read_forcefield.cpp: read_forcefield and ligand_stuff.cpp: read_ligand_data.cpp

*/

#include "pairwise_energy_calc.h"

void intrinsic_rotamer_energy(RESPARAM *resparam)
{
  BACKBONE bkbn;
  SIDECHAIN side;
  mini_pdbATOM *side_pdb, *dummyside_pdb,  CB_bkbn_pdb;
  mini_pdbATOM *firstAtom, *contactAtom, *twobondAtom, *dihedAtom;
  RESPARAM *resparam0;
  int q,i,pdb_ctr,pdb_ctr2, flag, step, doh,j,k;
  double r2; /* distance-squared */
  double r;
  double E_vdw, E_coulomb, E_vdw_14, E_coulomb_14, E_torsion;
  extern BACKBONE DUMMY_BKBN, PROLINE_BKBN;
  VDW_ENERGY vdwE;
  extern int INTRA_ROTAMER_FLAG_14, INTRA_ROTAMER_FLAG; 

  q=1;
    
  firstAtom = NULL; contactAtom = NULL; twobondAtom = NULL; dihedAtom = NULL;

  side_pdb = (mini_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(mini_pdbATOM));
  dummyside_pdb = (mini_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(mini_pdbATOM));
  side.atom = (micro_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(micro_pdbATOM));
  
  resparam0 = resparam; 
  ++resparam;
  doh=1;
  while(strcmp(resparam->residuetype,"zzz")!=0) /* loop over all residues */
    {
      if(abs(resparam->residuecode)<100)	/* pseudo-residues not considered; they don't have rotamers */
	{
	  
	  bkbn = DUMMY_BKBN; 
	  if(strcmp(resparam->residuetype,"PRO")==0)
	    bkbn = PROLINE_BKBN;
	    
	  CB_bkbn_pdb.coord = bkbn.CB;	    /* make CB atom */
	  CB_bkbn_pdb.atom_ptr = &resparam->atom[8];
	  CB_bkbn_pdb.seq_position=1;
	  CB_bkbn_pdb.sasa=0;

	 
	  for(i=1;i<=resparam->rotamerlib_ptr->number_of_rotamers;++i)
		resparam->rotamerlib_ptr->rotamer[i].coulombic = NULL;
	 
	  i=1;
	  while(i<=resparam->rotamerlib_ptr->number_of_rotamers)	/* loop over all rotamers */
	    {

	    if(resparam->rotamerlib_ptr->rotamer[i].coulombic == NULL) // if not NULL, this has already been done 
	      {	
		
		/* generate coordinates for this rotamer */

	      build_a_SideChain(&side, &bkbn,resparam,&q,resparam->rotamerlib_ptr->rotamer[i].chi);
	      resparam->rotamerlib_ptr->rotamer[i].born_factor = (double *)calloc((resparam->numberofAtoms+2), sizeof(double));
	      resparam->rotamerlib_ptr->rotamer[i].coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC));
	      resparam->rotamerlib_ptr->rotamer[i].coulombic->next = NULL;
	      resparam->rotamerlib_ptr->rotamer[i].first_coulombic = resparam->rotamerlib_ptr->rotamer[i].coulombic;

	      make_side_mini_pdbATOM(dummyside_pdb,&side);
	      pdb_ctr = 1;
	      side_pdb[pdb_ctr] = CB_bkbn_pdb; /* dummyside_pdb from make_side_mini_pdbATOM does not include CB, so add it in */
	      ++pdb_ctr;
	      while(dummyside_pdb[pdb_ctr-1].seq_position!=ENDFLAG)
	      {
	        side_pdb[pdb_ctr] = dummyside_pdb[pdb_ctr-1];
		++pdb_ctr;
	      }	
	      side_pdb[pdb_ctr].sasa = ENDFLAG;
	      side_pdb[pdb_ctr].seq_position = ENDFLAG;
	      side_pdb[pdb_ctr].atom_ptr=NULL; 

		/* zero the energies */
	      pdb_ctr = 1; /* E=0; */ 
	      ENERGY_0(resparam->rotamerlib_ptr->rotamer[i].energy);  VDW_ENERGY_0(vdwE);
	      E_vdw = 0; E_coulomb = 0; E_vdw_14 = 0; E_coulomb_14 = 0; E_torsion = 0; 

		/* start born factors w/ born_bonded */
	      while(side_pdb[pdb_ctr].seq_position!=ENDFLAG)
		{
		  resparam->rotamerlib_ptr->rotamer[i].born_factor[pdb_ctr] = 0;
		  if(side_pdb[pdb_ctr].atom_ptr->charge!=0)      
		    resparam->rotamerlib_ptr->rotamer[i].born_factor[pdb_ctr] = side_pdb[pdb_ctr].atom_ptr->born_bonded; 
		  ++pdb_ctr;
		}
	      resparam->rotamerlib_ptr->rotamer[i].born_factor[pdb_ctr]=ENDFLAG;
	      
		/* loop over all atom-pairs for this rotamer */
	      pdb_ctr = 1; 
	      while(side_pdb[pdb_ctr].seq_position!=ENDFLAG)
		{
		  if(side_pdb[pdb_ctr].atom_ptr->atom_ptr->atom_class != ENDFLAG)
		  {
		  pdb_ctr2 = pdb_ctr+1; 
		  while(side_pdb[pdb_ctr2].seq_position!=ENDFLAG)
		  {  
		      if(side_pdb[pdb_ctr2].atom_ptr->atom_ptr->atom_class != ENDFLAG)
		      {
		      r2 = Distance_sqrd(side_pdb[pdb_ctr].coord, side_pdb[pdb_ctr2].coord); 
		      if(r2>=TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD)
			{
			  r=sqrt(r2); /* E=0; */ flag=0;
	
				/* Determine connectivity */

			  /* atoms not directly bonded */
			  if((side_pdb[pdb_ctr].atom_ptr->other_info!=side_pdb[pdb_ctr2].atom_ptr->other_info) &&
			     (side_pdb[pdb_ctr].atom_ptr->contactAtomNumber != side_pdb[pdb_ctr2].atom_ptr->atomnumber) &&
			     (side_pdb[pdb_ctr2].atom_ptr->contactAtomNumber != side_pdb[pdb_ctr].atom_ptr->atomnumber) &&
			     (side_pdb[pdb_ctr].atom_ptr->ringContactAtom != side_pdb[pdb_ctr2].atom_ptr->atomnumber) &&
			     (side_pdb[pdb_ctr2].atom_ptr->ringContactAtom != side_pdb[pdb_ctr].atom_ptr->atomnumber))  
			    {	 
			      /* atoms not bonded to the same atom */
			      if((side_pdb[pdb_ctr].atom_ptr->contactAtomNumber != side_pdb[pdb_ctr2].atom_ptr->contactAtomNumber) &&
				 (side_pdb[pdb_ctr].atom_ptr->twobondAtomNumber != side_pdb[pdb_ctr2].atom_ptr->atomnumber) &&
				 (side_pdb[pdb_ctr2].atom_ptr->twobondAtomNumber != side_pdb[pdb_ctr].atom_ptr->atomnumber) &&
				 (side_pdb[pdb_ctr].atom_ptr->ringTwoBondAtom != side_pdb[pdb_ctr2].atom_ptr->atomnumber) &&
				 (side_pdb[pdb_ctr2].atom_ptr->ringTwoBondAtom != side_pdb[pdb_ctr].atom_ptr->atomnumber))
				{
				  if((side_pdb[pdb_ctr].atom_ptr->dihedAtomNumber == side_pdb[pdb_ctr2].atom_ptr->atomnumber) ||
				     (side_pdb[pdb_ctr2].atom_ptr->dihedAtomNumber == side_pdb[pdb_ctr].atom_ptr->atomnumber) ||
				     (side_pdb[pdb_ctr].atom_ptr->twobondAtomNumber == side_pdb[pdb_ctr2].atom_ptr->contactAtomNumber) ||
				     (side_pdb[pdb_ctr2].atom_ptr->twobondAtomNumber == side_pdb[pdb_ctr].atom_ptr->contactAtomNumber) ||
				     (side_pdb[pdb_ctr].atom_ptr->ringDihedAtom == side_pdb[pdb_ctr2].atom_ptr->atomnumber) ||
				     (side_pdb[pdb_ctr2].atom_ptr->ringDihedAtom == side_pdb[pdb_ctr].atom_ptr->atomnumber))
				    {    /* 1_4 interaction */
				      flag=INTRA_ROTAMER_FLAG_14;
				      if(pdb_ctr==1)
					flag=1;
				      if((side_pdb[pdb_ctr].atom_ptr->ringDihedAtom == side_pdb[pdb_ctr2].atom_ptr->atomnumber) || 
					 (side_pdb[pdb_ctr2].atom_ptr->ringDihedAtom == side_pdb[pdb_ctr].atom_ptr->atomnumber))
					{                                          /* through ring */
					  if(side_pdb[pdb_ctr].atom_ptr->ringDihedAtom == side_pdb[pdb_ctr2].atom_ptr->atomnumber)
					    {
					      firstAtom = &side_pdb[pdb_ctr];
					      dihedAtom = &side_pdb[pdb_ctr2];
					    }
					  else
					    {
					      firstAtom = &side_pdb[pdb_ctr2];
					      dihedAtom = &side_pdb[pdb_ctr];
					    }   
					  
					  if(firstAtom->atom_ptr->ringContactAtom<firstAtom->atom_ptr->atomnumber)
					    step=-1;
					  else
					    step=1;
					  contactAtom = firstAtom;
					  find_sidechain_atom(contactAtom,firstAtom->atom_ptr->ringContactAtom, step); 
					  
					  if(firstAtom->atom_ptr->ringTwoBondAtom<firstAtom->atom_ptr->atomnumber)
					    step=-1;
					  else
					    step=1;
					  twobondAtom = firstAtom;
					  find_sidechain_atom(twobondAtom, firstAtom->atom_ptr->ringTwoBondAtom, step);
					  
					}  
				      else      /* conventional 1_4 interaction */
					{	 
					  if((side_pdb[pdb_ctr].atom_ptr->dihedAtomNumber == side_pdb[pdb_ctr2].atom_ptr->atomnumber) ||
					     (side_pdb[pdb_ctr2].atom_ptr->dihedAtomNumber == side_pdb[pdb_ctr].atom_ptr->atomnumber))
					    {
					      if(side_pdb[pdb_ctr].atom_ptr->dihedAtomNumber == side_pdb[pdb_ctr2].atom_ptr->atomnumber)
						{
						  firstAtom = &side_pdb[pdb_ctr];
						  dihedAtom = &side_pdb[pdb_ctr2];
						}
					      else
						{
						  firstAtom = &side_pdb[pdb_ctr2];
						  dihedAtom = &side_pdb[pdb_ctr];
						} 
					      if(firstAtom->atom_ptr->contactAtomNumber<firstAtom->atom_ptr->atomnumber)
						step=-1;
					      else
						step=1;
					      contactAtom = firstAtom;
					      find_sidechain_atom(contactAtom,firstAtom->atom_ptr->contactAtomNumber, step); 
					      
					      if(firstAtom->atom_ptr->twobondAtomNumber<firstAtom->atom_ptr->atomnumber)
						step=-1;
					      else
						step=1;
					      twobondAtom = firstAtom;
					      find_sidechain_atom(twobondAtom, firstAtom->atom_ptr->twobondAtomNumber, step);		   
					    }
					  else
					    if((side_pdb[pdb_ctr].atom_ptr->contactAtomNumber == side_pdb[pdb_ctr2].atom_ptr->twobondAtomNumber) ||
					       (side_pdb[pdb_ctr2].atom_ptr->contactAtomNumber == side_pdb[pdb_ctr].atom_ptr->twobondAtomNumber))
					      {
						if(side_pdb[pdb_ctr].atom_ptr->contactAtomNumber == side_pdb[pdb_ctr2].atom_ptr->twobondAtomNumber)
						  {
						    firstAtom = &side_pdb[pdb_ctr2];
						    dihedAtom = &side_pdb[pdb_ctr];
						  }
						else
						  {
						    firstAtom = &side_pdb[pdb_ctr];
						    dihedAtom = &side_pdb[pdb_ctr2];
						  }
						
						if(firstAtom->atom_ptr->contactAtomNumber<firstAtom->atom_ptr->atomnumber)
						  step=-1;
						else
						  step=1;
						contactAtom = firstAtom;
						find_sidechain_atom(contactAtom,firstAtom->atom_ptr->contactAtomNumber, step); 
						
						if(firstAtom->atom_ptr->twobondAtomNumber<firstAtom->atom_ptr->atomnumber)
						  step=-1;
						else
						  step=1;
						twobondAtom = firstAtom;
						find_sidechain_atom(twobondAtom, firstAtom->atom_ptr->twobondAtomNumber, step);
					      }
					}
				    }	
				  else  /* atoms are greater than 3 bonds away */
				    {
				      flag=INTRA_ROTAMER_FLAG;	/* user-defined for calculating internal energies; default is to */
				    }
				}				     
			    } 

			  /* calculate >= 1_4 born factors */
				if(flag!=0)	/* >=1,4 connection; calculate desolvation */
				    {
				      if(side_pdb[pdb_ctr].atom_ptr->charge!=0)
				      {
				       if(strcmp(side_pdb[pdb_ctr2].atom_ptr->atom_ptr->atomtype, "HU")!=0)
					resparam->rotamerlib_ptr->rotamer[i].born_factor[pdb_ctr] += born_nonbond_calc(side_pdb[pdb_ctr].atom_ptr, side_pdb[pdb_ctr2].atom_ptr, r2); 
				      }
				      if(side_pdb[pdb_ctr2].atom_ptr->charge!=0)
				      {
				       if(strcmp(side_pdb[pdb_ctr].atom_ptr->atom_ptr->atomtype, "HU")!=0)
					resparam->rotamerlib_ptr->rotamer[i].born_factor[pdb_ctr2] += born_nonbond_calc(side_pdb[pdb_ctr2].atom_ptr, side_pdb[pdb_ctr].atom_ptr, r2);  
				      }
				    } 

				 /* store pair interactions in coulombic linked-list for born energy calcs;
						we don't have Born radii yet, since this is just an isolated rotamer.
					When the approximate born radii become available, go through this list and calculate
					the energies (in generate_lookup_table). */

				 if(flag<=3)    
				    if(side_pdb[pdb_ctr].atom_ptr->charge!=0) 
					if(side_pdb[pdb_ctr2].atom_ptr->charge!=0)  
					    {
						resparam->rotamerlib_ptr->rotamer[i].coulombic->index_i = pdb_ctr;
						resparam->rotamerlib_ptr->rotamer[i].coulombic->index_j = pdb_ctr2; 
						resparam->rotamerlib_ptr->rotamer[i].coulombic->r2 = r2;
						
						
						resparam->rotamerlib_ptr->rotamer[i].coulombic->flag = flag;
					    
						resparam->rotamerlib_ptr->rotamer[i].coulombic->next = (COULOMBIC *)malloc(sizeof(COULOMBIC));
						resparam->rotamerlib_ptr->rotamer[i].coulombic = resparam->rotamerlib_ptr->rotamer[i].coulombic->next;
						resparam->rotamerlib_ptr->rotamer[i].coulombic->r2 = ENDFLAG; /* mark end of COULOMBIC list */
						resparam->rotamerlib_ptr->rotamer[i].coulombic->next = NULL;
					    }   

 	  
			  if(flag!=0)
			    {
			      if(r2>=TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD)  
				{
				  if(flag==2)	/* vdw energy */
				    {
						add_atompair_vdw_energy(&vdwE, side_pdb[pdb_ctr].atom_ptr, side_pdb[pdb_ctr2].atom_ptr, r2, r);
						E_coulomb += coulomb_energy_q2_per_angstrom(r,side_pdb[pdb_ctr].atom_ptr,side_pdb[pdb_ctr2].atom_ptr);
				    }	 
				  if(flag==1)	/* 1,4 energy */
				    {	   
				      E_vdw_14 += vdw_energy_14(side_pdb[pdb_ctr].atom_ptr, side_pdb[pdb_ctr2].atom_ptr, r2);;
				      E_coulomb_14 += coulomb_energy_q2_per_angstrom(r,side_pdb[pdb_ctr].atom_ptr,side_pdb[pdb_ctr2].atom_ptr);
				      E_torsion += torsion_energy(firstAtom, contactAtom, twobondAtom, dihedAtom);
				   }	
				}
			    }
			}
		       }		      
		      ++pdb_ctr2;
		    }
		    }
		
		  ++pdb_ctr;
		}

	      resparam->rotamerlib_ptr->rotamer[i].coulombic = resparam->rotamerlib_ptr->rotamer[i].first_coulombic;
	      resparam->rotamerlib_ptr->rotamer[i].energy.E_vdw = sum_VDW_ENERGY(vdwE);
	      resparam->rotamerlib_ptr->rotamer[i].energy.E_coulomb = COULOMB_CONST_OVER_INTERNAL_DIELECTRIC*E_coulomb;
	      resparam->rotamerlib_ptr->rotamer[i].energy.E_1_4 = WEIGHT_1_4*(E_torsion + NONBOND_FACTOR_1_4*(E_vdw_14 + COULOMB_CONST*E_coulomb_14));
	    
	      	
		if(resparam->ligand_flag == 1)
		{
		if(resparam->num_rotatable_bonds == 0)	// for a rigid ligand, internal energies are all the same
		{
			for(i=2;i<=resparam->rotamerlib_ptr->number_of_rotamers;++i)
			{
				resparam->rotamerlib_ptr->rotamer[i].first_coulombic = resparam->rotamerlib_ptr->rotamer[1].first_coulombic;
				resparam->rotamerlib_ptr->rotamer[i].coulombic = resparam->rotamerlib_ptr->rotamer[i].first_coulombic;
				resparam->rotamerlib_ptr->rotamer[i].energy = resparam->rotamerlib_ptr->rotamer[1].energy;
				resparam->rotamerlib_ptr->rotamer[i].born_factor = resparam->rotamerlib_ptr->rotamer[1].born_factor;
			}

			i = resparam->rotamerlib_ptr->number_of_rotamers + 10;
		}
		else // flexible ligand; ligamers w/ same rot bonds are all the same 
		{
			// find other ligamers w/ the same rotamer and set up their pointers to this ligamer
			
			for(j=(i+1);j<=resparam->rotamerlib_ptr->number_of_rotamers;++j)
			{
					k=13;
					while(k<=resparam->rotamerlib_ptr->numChi)
					{
						if(resparam->rotamerlib_ptr->rotamer[i].chi[k] != resparam->rotamerlib_ptr->rotamer[j].chi[k])
							k = resparam->rotamerlib_ptr->numChi + 1000;

						++k;
					}
					--k;
					
					if(k == resparam->rotamerlib_ptr->numChi)  // ligamer j has the same rotamer as ligamer i
					{
				resparam->rotamerlib_ptr->rotamer[j].first_coulombic = resparam->rotamerlib_ptr->rotamer[i].first_coulombic;
				resparam->rotamerlib_ptr->rotamer[j].coulombic = resparam->rotamerlib_ptr->rotamer[j].first_coulombic;
				resparam->rotamerlib_ptr->rotamer[j].energy = resparam->rotamerlib_ptr->rotamer[i].energy;
				resparam->rotamerlib_ptr->rotamer[j].born_factor = resparam->rotamerlib_ptr->rotamer[i].born_factor;
					}

			}			
		}
		}

		} // end coulombic = NULL
		
		++i;
	    } // end i
	} // end abs(resparam->residuecode)<100
      
      ++resparam; ++doh;
    }

  
  resparam = resparam0; 
  
  free_memory(side_pdb);
  free_memory(dummyside_pdb);
  free_memory(side.atom);

  

}



/*
  Given the four atoms, calculate the dihedral-angle-dependent energy component of 1,4 interactions

  firstAtom
           \
	    contactAtom--------twobondAtom
					       \
					       dihedAtom
 */

double torsion_energy(mini_pdbATOM *firstAtom, mini_pdbATOM *contactAtom, mini_pdbATOM *twobondAtom, mini_pdbATOM *dihedAtom)	
{
    double phi;
    
    if(firstAtom->atom_ptr->atom_ptr->torsion !=NULL)	
    {
    phi = dihedral(&firstAtom->coord,&contactAtom->coord,&twobondAtom->coord,&dihedAtom->coord);
   
	firstAtom->atom_ptr->atom_ptr->torsion = firstAtom->atom_ptr->atom_ptr->first_torsion;
	
		/* go through linked list of torsions that contain this atomtype until the correct one is found */
	    while(firstAtom->atom_ptr->atom_ptr->torsion->next!=NULL)  
	    {
		if(firstAtom->atom_ptr->atom_ptr->torsion->torsion_class[2] == contactAtom->atom_ptr->atom_ptr->torsion_class)
		{
		    if(firstAtom->atom_ptr->atom_ptr->torsion->torsion_class[3] == twobondAtom->atom_ptr->atom_ptr->torsion_class)
		    {
			if(firstAtom->atom_ptr->atom_ptr->torsion->torsion_class[4] == dihedAtom->atom_ptr->atom_ptr->torsion_class)
			{    	
				   return( ((firstAtom->atom_ptr->atom_ptr->torsion->magnitude[1]*(1 + cosd(phi)) + 
					firstAtom->atom_ptr->atom_ptr->torsion->magnitude[2]*(1 - cosd(2.0*phi)) + 
					firstAtom->atom_ptr->atom_ptr->torsion->magnitude[3]*(1 + cosd(3.0*phi)))) ); 
					
			}
		    }		
		}	
		firstAtom->atom_ptr->atom_ptr->torsion = firstAtom->atom_ptr->atom_ptr->torsion->next;
	    }
    }
    else	/* no angle-dependent torsion energies defined for firstAtom's atomtype */
    {
	if(strcmp(firstAtom->atom_ptr->atom_ptr->torsiontype, "U")!=0)
	{
	fprintf(stderr,"ERROR torsion not found:\t");
	fprintf(stderr,"%d\t%s\t%s\t%s\t%s\t%d\n",  firstAtom->seq_position, 
				    firstAtom->atom_ptr->atom_ptr->atomtype, 
				    contactAtom->atom_ptr->atom_ptr->atomtype,
				    twobondAtom->atom_ptr->atom_ptr->atomtype,
				    dihedAtom->atom_ptr->atom_ptr->atomtype,
				    dihedAtom->seq_position);
	exit(1);
	}			    
   }

 /*  firstAtom's atomtype has angle-dependent torsion energies defined, but not with the other atoms inputed here */	   
   return(0.0); 
}

/* free a COULOMBIC linked-list; standard recursive linked-list freeing algorithm */
void freeCOULOMBIC(COULOMBIC *coulombic)
{
    COULOMBIC *temp;
    
    while(coulombic!=NULL)  
    { 
	temp=coulombic->next; 
	free_memory(coulombic);     
	coulombic=temp;      
    }
}

/* calculates the desolvation components for a pair of non-bonded (>=1,4) atoms for Born radii calcs
	GB implementation of Dominy and Brooks J Phys Chem B 103:3765-3773 (1999) */
double born_nonbond_calc(ATOMRESPARAM *atom_i, ATOMRESPARAM *atom_j, double r2)	 
{    
    double ccf, ratio;
    extern double ONE_OVER_BORN_P5, PI_TIMES_BORN_P5;
    extern double **UNITED_RADII_SUM_SQUARED;

    ratio = r2/UNITED_RADII_SUM_SQUARED[atom_i->atom_ptr->atom_class][atom_j->atom_ptr->atom_class];   
    
    if(ratio > ONE_OVER_BORN_P5)    
	return(atom_j->BORN_P4_volume/(r2*r2));	    
    else
    {
	ccf = 0.5*(1.0 - cos(ratio*PI_TIMES_BORN_P5));
	return(atom_j->BORN_P4_volume*ccf*ccf/(r2*r2));
    }
}
 

/* 
   This function calculates all energies for the structure in pdb; returned as ENERGY.
   pdb must have atom_ptrs assigned and be terminated with atom_ptr=NULL, seq_position=ENDFLAG.
   electr_calc_flag=0: no born calculations
   electr_calc_flag=1: calculate born radii and electrostatics.
   electr_calc_flag=2: calculate electrostatics w/ born radii already in pdb[i].born_radius.
   electr_calc_flag=3: store born desolvation sums in pdb[i].born_radius, but don't calculate born/pol energy
   sasa_flag=1: calculate surface-area based energies.	
	surface area components also stored in extern SASA_SUM SASA_SUM_TEMP; 
	this can be accessed by the calling function
*/
   
ENERGY energy_calc(mini_pdbATOM *pdb, int electr_calc_flag, int input_sasa_flag) 
{
  double r2; /* squared distance between two atoms */
  double EE, r;
  ENERGY energy;
  double E_coulomb, E_coulomb_14, E_vdw_14, E_torsion, E_hbond;
  double  alpha, f;
  double xmax, ymax, zmax, xmin, ymin, zmin;
  int flag2,  i_ctr, j_ctr;
  int i_heavy, j_heavy;
  mini_pdbATOM *firstAtom, *dihedAtom, *twobondAtom, *contactAtom;
  COULOMBIC *first_coulombic, *coulombic;
  double dx, dy, dz;
  int i, j, free_vector_ptr_flag, sasa_flag;
  CARTESIAN_VECTOR *i_j_vector, *dummy_vector;
  mini_pdbATOM *micro_ptr[2];
  SASA_SUM sasa_sum;
  extern CARTESIAN_VECTOR ***VECTOR_PTR;
  VDW_ENERGY vdwE;
  extern SASA_SUM SASA_SUM_TEMP;
  extern double **CHARGE_PRODUCT;

  VDW_ENERGY_0(vdwE);

  firstAtom = NULL; contactAtom = NULL; twobondAtom = NULL; dihedAtom = NULL;
  first_coulombic = NULL; coulombic = NULL;  i_j_vector = NULL; dummy_vector = NULL;

	/* initialize COULOMBIC linked-list */
  if(electr_calc_flag != 0 && electr_calc_flag != 3)
  {
    coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC));
    first_coulombic = coulombic;   
    coulombic->next = NULL;  
  }
  
  ENERGY_0(energy);
  E_coulomb=0;  E_coulomb_14=0; E_vdw_14=0; E_torsion=0;  E_hbond = 0;
  
	

  if(input_sasa_flag == 1)		/* prevents crashes due to CARTESIAN_VECTOR matrix being too big */
  {
  	if(MAX_ATOMS <= REAL_MAX_ATOMS)
		sasa_flag =1;
  	else
		sasa_flag =2;
  }
  else
	sasa_flag = input_sasa_flag;

/* matrix for SASA calculations; don't use if the structure is too big */
  free_vector_ptr_flag = 0;
  if(sasa_flag==1)
  {			    /* if this is NULL, then it means that this matrix should be created and destroyed internally */
    if(VECTOR_PTR==NULL)    /* for applications that require many calls to this function, this matrix should be allocated elsewhere */
    {			    /* and kept for the duration of the need; this matrix is the largest object in terms of memory */
	free_vector_ptr_flag=1;	
		
	VECTOR_PTR = (CARTESIAN_VECTOR ***)calloc(MAX_ATOMS, sizeof(CARTESIAN_VECTOR **));
	for(i=0;i<MAX_ATOMS;++i)
	    VECTOR_PTR[i] = (CARTESIAN_VECTOR **)calloc(MAX_ATOMS, sizeof(CARTESIAN_VECTOR *));
    }
  
    i_j_vector = (CARTESIAN_VECTOR *)malloc(sizeof(CARTESIAN_VECTOR));
    dummy_vector = i_j_vector;
    i_j_vector->next = NULL;
    micro_ptr[0] = &pdb[1];
    micro_ptr[1] = NULL;
  }  
  
    i_ctr = 1; 
    while(pdb[i_ctr].seq_position!=ENDFLAG)
    {
	pdb[i_ctr].hbond_satisfied = 1;

	if(pdb[i_ctr].atom_ptr->donorflag[0]!=' ')
	{
			pdb[i_ctr].hbond_satisfied = 0;	/* will check for satisfaction later */
		
	}
		/* initialize born_radius of charged atoms w/ born_bonded if electr_calc_flag==1 or 3 */
		/* if 2, use born radii already in pdb[i_ctr].born_radius */
	if(electr_calc_flag==1 || electr_calc_flag==3)
		if(pdb[i_ctr].atom_ptr->charge!=0)
			if(pdb[i_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG)
				pdb[i_ctr].born_radius = pdb[i_ctr].atom_ptr->born_bonded; 

	++i_ctr;
    }
 
  
  i_ctr = 1; i_heavy=0;
  while(pdb[i_ctr].seq_position!=ENDFLAG)   /* i loop */
    {
    if(pdb[i_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG)  /* no U */
    {
	
      if(sasa_flag==1 && pdb[i_ctr].atom_ptr->atom_ptr->sasa_radius>0)
	{
	  ++i_heavy;	/* initialize box for sasa calculations */
	  pdb[i_ctr].box = ((int)(pdb[i_ctr].coord.x/SASA_DISTANCE_CUTOFF)+500)*1000000	
	    +((int)(pdb[i_ctr].coord.y/SASA_DISTANCE_CUTOFF)+500)*1000
	    +(int)(pdb[i_ctr].coord.z/SASA_DISTANCE_CUTOFF)+500
	    -(pdb[i_ctr].coord.x<0)*1000000
	    -(pdb[i_ctr].coord.y<0)*1000
	    -(pdb[i_ctr].coord.z<0);
	}	
      
	  xmax=pdb[i_ctr].coord.x+FORCEFIELD_DISTANCE_CUTOFF;
	  xmin=pdb[i_ctr].coord.x-FORCEFIELD_DISTANCE_CUTOFF;
	  ymax=pdb[i_ctr].coord.y+FORCEFIELD_DISTANCE_CUTOFF;
	  ymin=pdb[i_ctr].coord.y-FORCEFIELD_DISTANCE_CUTOFF;
	  zmax=pdb[i_ctr].coord.z+FORCEFIELD_DISTANCE_CUTOFF;
	  zmin=pdb[i_ctr].coord.z-FORCEFIELD_DISTANCE_CUTOFF;
	  
      j_ctr= i_ctr+1; j_heavy = i_heavy;
      while(pdb[j_ctr].seq_position!=ENDFLAG)  /* j loop */
	{
	if(pdb[j_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG)
	{
	  if(pdb[j_ctr].atom_ptr->atom_ptr->sasa_radius>0)
	    ++j_heavy;
	  /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF
	     centered at pdb[j_ctr].coord */
	  if(pdb[j_ctr].coord.x<=xmax)
	    if(pdb[j_ctr].coord.x>=xmin)
	      if(pdb[j_ctr].coord.y<=ymax)
		if(pdb[j_ctr].coord.y>=ymin)
		  if(pdb[j_ctr].coord.z<=zmax)
		    if(pdb[j_ctr].coord.z>=zmin)
		      {
			dx = pdb[j_ctr].coord.x - pdb[i_ctr].coord.x;
			dy = pdb[j_ctr].coord.y - pdb[i_ctr].coord.y;
			dz = pdb[j_ctr].coord.z - pdb[i_ctr].coord.z;
			r2 = dx*dx + dy*dy + dz*dz;  
			if(r2>=TINY &&  r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD)  /* FORCEFIELD_DISTANCE_CUTOFF cutoff */  
			  {
			    flag2=0;  
			    r=sqrt(r2);
			    if(sasa_flag==1)
			      {	 
				if(pdb[i_ctr].atom_ptr->atom_ptr->sasa_radius>0)
				  if(pdb[j_ctr].atom_ptr->atom_ptr->sasa_radius>0)
				    {

				/* if these atoms bury each other, measure the interbond vector for sasa calcs */
				/* triangle inequality; pointed out by Mark Voorhies */
			  if(r2<=SASA_RADII_SQRD_MATRIX[pdb[i_ctr].atom_ptr->atom_ptr->atom_class][pdb[j_ctr].atom_ptr->atom_ptr->atom_class])
				if(r + pdb[i_ctr].atom_ptr->atom_ptr->sasa_radius >= pdb[j_ctr].atom_ptr->atom_ptr->sasa_radius)
				if(r + pdb[j_ctr].atom_ptr->atom_ptr->sasa_radius >= pdb[i_ctr].atom_ptr->atom_ptr->sasa_radius)
					{
					  VECTOR_PTR[i_heavy][j_heavy] = i_j_vector;
					  VECTOR_PTR[j_heavy][i_heavy] = i_j_vector;
					  i_j_vector->i = i_heavy;
					  i_j_vector->j = j_heavy;
					  i_j_vector->r2 = r2;
					  i_j_vector->r = r;
					  dx/=r; dy/=r; dz/=r;
					  i_j_vector->i_theta = acos(dz);
					  i_j_vector->i_phi = atan2(dx,dy);
					  
					  dx = -dx; dy= -dy; dz= -dz;
					  i_j_vector->j_theta = acos(dz);
					  i_j_vector->j_phi = atan2(dx,dy);
					  
					  i_j_vector->next =  (CARTESIAN_VECTOR *)malloc(sizeof(CARTESIAN_VECTOR));
					  i_j_vector = i_j_vector->next;	
					  i_j_vector->next =  NULL;
					}	
				    }
			      }
			      
					/* get the connectivity between the atoms; if 1,4, also get pointers to the dihedral atoms */
				flag2 = get_atom_connectivity(r2, &pdb[i_ctr], &pdb[j_ctr], 
												&firstAtom, &contactAtom, &twobondAtom, &dihedAtom);

			/* calculate desolvation for born radii calcs  */
			     if(pdb[i_ctr].atom_ptr->charge!=0) /*  i charged  */
			      {   
			          if(flag2!=0) 
				    if(electr_calc_flag == 1 || electr_calc_flag == 3)    
				     if(strcmp(pdb[j_ctr].atom_ptr->atom_ptr->atomtype, "HU")!=0)  /* HU does not influence the born radius */
					pdb[i_ctr].born_radius += born_nonbond_calc(pdb[i_ctr].atom_ptr, pdb[j_ctr].atom_ptr, r2); 
				     
				if(pdb[j_ctr].atom_ptr->charge!=0)  /* i and j charged */
				  {
				    if(flag2!=0) 
				      if(electr_calc_flag == 1 || electr_calc_flag == 3)
					if(strcmp(pdb[i_ctr].atom_ptr->atom_ptr->atomtype, "HU")!=0)  /* HU does not influence the born radius */
					  pdb[j_ctr].born_radius += born_nonbond_calc(pdb[j_ctr].atom_ptr, pdb[i_ctr].atom_ptr, r2); 
				    
					
				    if(flag2<=3)	/* save distances for born calcs later */
				       {
					    if(electr_calc_flag != 0 && electr_calc_flag != 3) 
						{
						coulombic->index_i = i_ctr;
						coulombic->index_j = j_ctr; 
						coulombic->r2 = r2;
						coulombic->flag = flag2;
					    
						coulombic->next = (COULOMBIC *)malloc(sizeof(COULOMBIC));  /* allocate memory for next element */
						coulombic = coulombic->next;  
						coulombic->next = NULL; /* last element unless replaced */
					    

						}
					  } 
				    }  		     
			      }  	
			    else
			     if(pdb[j_ctr].atom_ptr->charge!=0) /*  only j charged  */
				{    
				  if(flag2!=0) 
				    if(electr_calc_flag == 1 || electr_calc_flag == 3)
				      {
				        if(strcmp(pdb[i_ctr].atom_ptr->atom_ptr->atomtype, "HU")!=0)  /* HU does not influence the born radius */
					    pdb[j_ctr].born_radius += born_nonbond_calc(pdb[j_ctr].atom_ptr, pdb[i_ctr].atom_ptr, r2);
				      }
			   	}
			   
			     if(flag2!=0) 
			    {


				  if(flag2==2)    /*   >1,4 interaction  */
				    {
						E_hbond = hbond_energy(&pdb[i_ctr], &pdb[j_ctr], r, r2);
						if(E_hbond != 0)
						{
							energy.E_hbond += E_hbond;
							++energy.num_hbond;
							if(pdb[i_ctr].atom_ptr->other_info<0 && pdb[j_ctr].atom_ptr->other_info<0)
								++energy.num_bb_hbond;
							else if(pdb[i_ctr].atom_ptr->other_info>0 && pdb[j_ctr].atom_ptr->other_info>0)
								++energy.num_ss_hbond;
							else
								++energy.num_sb_hbond;
						}
						
						E_coulomb += coulomb_energy_q2_per_angstrom(r, pdb[i_ctr].atom_ptr, pdb[j_ctr].atom_ptr);
						
						add_atompair_vdw_energy(&vdwE, pdb[i_ctr].atom_ptr, pdb[j_ctr].atom_ptr, r2, r);

					/*
						if(DEBUG == 1)
						{

							VDW_ENERGY dummy_vdwE;
							VDW_ENERGY_0(dummy_vdwE);
							
							add_atompair_vdw_energy(&dummy_vdwE, pdb[i_ctr].atom_ptr, pdb[j_ctr].atom_ptr, r2, r);

							printf("%d\t%s\t%s\t%d\t%s\t%s\t%lf\t%lf\n",
								pdb[i_ctr].seq_position,pdb[i_ctr].atom_ptr->atomname,pdb[i_ctr].atom_ptr->atomtype,
								pdb[j_ctr].seq_position,pdb[j_ctr].atom_ptr->atomname,pdb[j_ctr].atom_ptr->atomtype,
								r,sum_VDW_ENERGY(dummy_vdwE));
						}
					*/	


						

				    } 
				   if(flag2==1)       /* 1,4 */
				    {
					
					  E_vdw_14 += vdw_energy_14(pdb[i_ctr].atom_ptr, pdb[j_ctr].atom_ptr, r2);
				      E_coulomb_14 += coulomb_energy_q2_per_angstrom(r, pdb[i_ctr].atom_ptr, pdb[j_ctr].atom_ptr);
				      E_torsion += torsion_energy(firstAtom, contactAtom, twobondAtom, dihedAtom); 
				     }   
			    }   
			}
		   }
	}    
	++j_ctr;
	}
    }

	energy.E_specificity += 
			hbonding_status_dependent_specificity_energy(pdb[i_ctr].hbond_satisfied, pdb[i_ctr].atom_ptr);
	

    ++i_ctr;
    }
   

	/* calculate SASA */
  if(sasa_flag != 0)
  {
  	if(sasa_flag==1)
    	{
    		sasa_sum = precalcSASAcalc(micro_ptr);
    	
    		i_j_vector = dummy_vector;
    		freeCARTESIAN_VECTOR(i_j_vector);
  
    		for(i=1;i<=i_heavy;++i)
			for(j=1;j<=i_heavy;++j)
			{
	    			if(VECTOR_PTR[i][j]!=NULL)
	    			{
					VECTOR_PTR[i][j]=NULL;
					VECTOR_PTR[j][i]=NULL;
	    			}    
			}
	
    		if(free_vector_ptr_flag==1)
    		{   
			for(i=0;i<MAX_ATOMS;++i)
	    			free_memory(VECTOR_PTR[i]);
			free_memory(VECTOR_PTR);  
			VECTOR_PTR = NULL;
    		}
   	}
   	else if(sasa_flag ==2)
  	{
		sasa_sum = SASAcalc(pdb);
  	} 
    
    	energy.E_sasa = GENERAL_ASP*sasa_sum.sasa_total + HYDROPHOBIC_ASP*(sasa_sum.sasa_total - (sasa_sum.N + sasa_sum.O));	
   
    	SASA_SUM_TEMP = sasa_sum; /* so the calling function can see the surface area components */
  }	


    /* calculate born energies */
  energy.E_born = 0; energy.E_pol=0;
  if(electr_calc_flag!=0 && electr_calc_flag!=3)
    {
	/* self energies for each atom, using the desolvation sum stored in pdb[i_ctr].born_radius 
		calculate born radii here */

	energy.E_specificity=0;

	if(sasa_flag != 0)
	{
		exposed_hbond_groups(pdb); 
	}

      i_ctr = 1; energy.num_unsatisfied_hbond=0; 
	energy.num_side_unsatisfied_hbond=0; energy.num_bkbn_unsatisfied_hbond=0; 
      while(pdb[i_ctr].seq_position!=ENDFLAG)  
	{
	  energy.E_specificity += 
			hbonding_status_dependent_specificity_energy(pdb[i_ctr].hbond_satisfied, pdb[i_ctr].atom_ptr);


	  if(pdb[i_ctr].atom_ptr->charge!=0)
	  if(pdb[i_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG)
	    {   
		  if(electr_calc_flag == 1)    
		  {
		    if(pdb[i_ctr].born_radius >= 0)
			pdb[i_ctr].born_radius = -1.0;
		    pdb[i_ctr].born_radius = -HALF_COULOMB_CONST/pdb[i_ctr].born_radius;
		  } 
		EE = pdb[i_ctr].atom_ptr->charge2/pdb[i_ctr].born_radius;
		
		energy.E_born += BORN_CONST(pdb[i_ctr].born_radius)*EE; 

		
		
			if(pdb[i_ctr].hbond_satisfied == 0)
			{
				++energy.num_unsatisfied_hbond;
		 		if(pdb[i_ctr].atom_ptr->other_info<0)
					++energy.num_bkbn_unsatisfied_hbond;
				else
					++energy.num_side_unsatisfied_hbond;
			}

	    }	  
	  ++i_ctr;
	}


 	// printf("%d\t",energy.num_unsatisfied_hbond);
	     
	/* born cross-polarization energies; use the coulombic linked-list with the born radii calculated above */
      coulombic = first_coulombic;
      while(coulombic->next!=NULL)
	{
	  alpha = pdb[coulombic->index_i].born_radius*pdb[coulombic->index_j].born_radius;
	   
	  f = f_of_r2_alpha(coulombic->r2, alpha);
	  
	  EE = CHARGE_PRODUCT[pdb[coulombic->index_i].atom_ptr->charge_class][pdb[coulombic->index_j].atom_ptr->charge_class]/f; 
	  
	/* prevent favorable born pair energies from overwhelming lennard-jones for close, nonbonded atoms */
	  if(coulombic->flag == 2)
	  {
	      if(EE > POL_ATTRACTION_CAP)	 
		EE = POL_ATTRACTION_CAP;
	  }

	  energy.E_pol += BORN_CONST(f)*EE;
	
	  coulombic = coulombic->next;      
	}
    
    coulombic = first_coulombic;
    freeCOULOMBIC(coulombic);

    coulombic = NULL; 
    first_coulombic = NULL;  
  }

	/* add up energies, multiply with constants to get the correct units, etc */  

   energy.E_coulomb = COULOMB_CONST_OVER_INTERNAL_DIELECTRIC*E_coulomb;

   
	energy.E_pol = COULOMB_CONST*energy.E_pol;   
	energy.E_born = HALF_COULOMB_CONST*energy.E_born;   
   
	   
  energy.E_1_4 = WEIGHT_1_4*(NONBOND_FACTOR_1_4*(E_vdw_14 + COULOMB_CONST*E_coulomb_14) + E_torsion);

  energy.E_vdw = sum_VDW_ENERGY(vdwE);	
  energy.clashes = vdwE.clashes; 

  energy.E_hbond = HBOND_FUNCTION_FLAG*energy.E_hbond;

 return energy;
  
}

/*	determines connectivity between two atoms pdb_i and pdb_j in an array of mini_pdbATOM in the calling function 
	returns 0 for <=1,3; 1 for 1,4; 2 for >1,4 ; for 1,4, returns non-NULL values for dihedral atoms */

int get_atom_connectivity(double r2, mini_pdbATOM *pdb_i, mini_pdbATOM *pdb_j, 
	mini_pdbATOM **firstAtom, mini_pdbATOM **contactAtom, mini_pdbATOM **twobondAtom, mini_pdbATOM **dihedAtom)
{
	int connection_flag;
	int step, stop;
	mini_pdbATOM *dummybkbn, *dummybkbn2, *dummy_side, *dummy_bkbn;
	extern int INTRA_ROTAMER_FLAG_14, INTRA_ROTAMER_FLAG; 

	*firstAtom=NULL;
	*contactAtom=NULL;
	*twobondAtom=NULL;
	*dihedAtom=NULL;

	connection_flag=0;

			/* backbone-backbone */
			    if(pdb_i->atom_ptr->other_info<0 && pdb_j->atom_ptr->other_info<0)
			      {
				if(pdb_i->seq_position==pdb_j->seq_position)   /* same residue */
				  {
				    if(strcmp(pdb_i->atom_ptr->atom_ptr->atomtype,"HN")==0 || 
				       strcmp(pdb_j->atom_ptr->atom_ptr->atomtype,"HN")==0)
				      {
					if(strcmp(pdb_i->atom_ptr->atom_ptr->atomtype,"HN")==0)  /* i_atom is H */
					  {
					    dummybkbn=pdb_i; dummybkbn2=pdb_j;
					  }
					else      /* j_atom is H */
					  {					    
					    dummybkbn2=pdb_i; dummybkbn=pdb_j;
					  }
					if(strstr(dummybkbn2->atom_ptr->atomname,"HA")!=0 ||
					   strcmp(dummybkbn2->atom_ptr->atomname,"CB")==0 ||
					   strcmp(dummybkbn2->atom_ptr->atomname,"C")==0)
					  {
					    connection_flag=1;  /* 1,4 */
					    *firstAtom = dummybkbn;
					    *dihedAtom = dummybkbn2;
					    
					    *contactAtom = *firstAtom;
					    find_backbone_atom((*contactAtom),"N", -1); 
					     
					    *twobondAtom = *firstAtom;
					    find_backbone_atom((*twobondAtom),"CA", 1); 
					    
					  }
					else
					  if(strstr(dummybkbn2->atom_ptr->atomname,"O")!=0)
					    connection_flag=2; /* >1,4 */
				      }   
				    else
				    {
				      if(strstr(pdb_i->atom_ptr->atomname,"O")!=0 || 
				         strstr(pdb_j->atom_ptr->atomname,"O")!=0)
					{
					  if(strstr(pdb_i->atom_ptr->atomname,"O")!=0)   /* i_atom is O */
					    {
					      dummybkbn=pdb_i; dummybkbn2=pdb_j;
					    }
					  else      /* j_atom is O */
					    {					    
					      dummybkbn2=pdb_i; dummybkbn=pdb_j;
					    }
					  if(strcmp(dummybkbn2->atom_ptr->atomname,"N")==0 || 
					     strstr(dummybkbn2->atom_ptr->atomname,"HA")!=0 ||
					     strcmp(dummybkbn2->atom_ptr->atomname,"CB")==0)
					    {
					      connection_flag=1;  /* 1,4 */
					      *firstAtom = dummybkbn;
					      *dihedAtom = dummybkbn2;
					      
					      *contactAtom = *firstAtom;
					      find_backbone_atom((*contactAtom), "C", -1);
					      
					      *twobondAtom = *contactAtom;
					      find_backbone_atom((*twobondAtom), "CA", -1);
					      	
					    }	
					  else
					    if(strcmp(dummybkbn2->atom_ptr->atom_ptr->atomtype,"HN")==0)
						connection_flag=2; 	
					}
				     else
					connection_flag=0;
				    }									  
				  }	  
				else
				  {
				    if(abs(pdb_i->seq_position - pdb_j->seq_position)==1) /* neighboring residue */
				      {
				        if(pdb_i->seq_position > pdb_j->seq_position)
					  {
					    dummybkbn=pdb_j; dummybkbn2=pdb_i;
					  }
				        else
					  {
					    dummybkbn=pdb_i; dummybkbn2=pdb_j;
					  }
					
					if(strcmp(dummybkbn->atom_ptr->atomname,"N")==0 || 
					   strcmp(dummybkbn->atom_ptr->atomname,"CB")==0 || 
					   strstr(dummybkbn->atom_ptr->atomname,"HA")!=0)
					  {
					    if(strcmp(dummybkbn2->atom_ptr->atomname,"N")==0)
					      {
						connection_flag=1;
						*firstAtom = dummybkbn;
						*dihedAtom = dummybkbn2;
						
						*contactAtom = dummybkbn;
						find_backbone_atom((*contactAtom), "CA", 1);
						
						*twobondAtom = *contactAtom;
						find_backbone_atom((*twobondAtom), "C", 1);
	
					      }  
					    else
					      connection_flag=2;
					  }
					else
					  if(strcmp(dummybkbn->atom_ptr->atom_ptr->atomtype,"HN")==0 ||
						 strcmp(dummybkbn->atom_ptr->atom_ptr->atomtype,"CD")==0 ||
						 strcmp(dummybkbn->atom_ptr->atom_ptr->atomtype,"CG")==0) 
					    {
					      connection_flag=2;
					    }
					  else
					    if(strcmp(dummybkbn->atom_ptr->atomname,"C")==0)
					      {
						if(strcmp(dummybkbn2->atom_ptr->atomname,"C")==0 || 
						   strcmp(dummybkbn2->atom_ptr->atomname,"CB")==0 ||
						   strstr(dummybkbn2->atom_ptr->atomname,"HA")!=0)
						  {
						    connection_flag=1;
						    *firstAtom = dummybkbn2;
						    *dihedAtom = dummybkbn;
						    
						    *contactAtom = dummybkbn2;
						    find_backbone_atom((*contactAtom), "CA", -1);
						    
						    *twobondAtom = *contactAtom;
						    find_backbone_atom((*twobondAtom), "N", -1);
						   
						  }
						else
						  if(strstr(dummybkbn2->atom_ptr->atomname,"O")!=0)
						    connection_flag=2;
						  else if(strcmp(dummybkbn2->atom_ptr->atomname,"CG")==0) /* PRO */
							{
								connection_flag=1;
								*firstAtom = dummybkbn2;
						    	*dihedAtom = dummybkbn;

								*contactAtom = dummybkbn2;
						    	find_backbone_atom((*contactAtom), "CD", 1);

								*twobondAtom = dummybkbn2;
						    	find_backbone_atom((*twobondAtom), "N", -1);

							}
					      }       
					    else
					      if(strcmp(dummybkbn->atom_ptr->atomname,"CA")==0 || 
					         strstr(dummybkbn->atom_ptr->atomname,"O")!=0)
						{
						  if(strcmp(dummybkbn2->atom_ptr->atomname,"CA")==0 || 
						     strcmp(dummybkbn2->atom_ptr->atom_ptr->atomtype,"HN")==0 ||
							 strcmp(dummybkbn2->atom_ptr->atom_ptr->atomtype,"CD")==0) /* PRO */
						    {
						      connection_flag=1;
						      *firstAtom = dummybkbn2;
						      *dihedAtom = dummybkbn;
						      
						      *contactAtom = dummybkbn2;
						      find_backbone_atom((*contactAtom), "N", -1);
						      
						      *twobondAtom = dummybkbn;
						      find_backbone_atom((*twobondAtom), "C", 1);
						      
						    }
						  else
						    if(strcmp(dummybkbn2->atom_ptr->atomname,"N")!=0)
						      connection_flag=2;   
						}
						else
						    connection_flag=2;    	   
				      }   
				    else                /* non-neighboring residue */
				      connection_flag=2;
				      	 
				  }
			      }
			    else   /* sidechain-sidechain or sidechain-bkbn */
			      {
				if(pdb_i->seq_position!=pdb_j->seq_position)  /* not within same residue */
				{ 
				     connection_flag=2;
				
					/* don't score disulfide-bonded sulfurs */
					if(r2<=MAX_DISULFIDE_BOND_LENGTH_SQRD)
					{
					  if(strcmp(pdb_i->atom_ptr->atom_ptr->atomtype,"SH")==0)
						if(strcmp(pdb_j->atom_ptr->atom_ptr->atomtype,"SH")==0)
							connection_flag=0;
					}
					
				 }   
				else   /* within same residue */
				  {
				    /* intra-sidechain */
				     if(pdb_i->atom_ptr->other_info>0 && pdb_j->atom_ptr->other_info>0)
					{
				 	 if(pdb_i->atom_ptr->other_info != pdb_j->atom_ptr->other_info) /* filter out within same aromatic */
					  {
					  /* atoms not directly bonded */
					  if((pdb_i->atom_ptr->contactAtomNumber != pdb_j->atom_ptr->atomnumber) &&
					     (pdb_j->atom_ptr->contactAtomNumber != pdb_i->atom_ptr->atomnumber) &&
					     (pdb_i->atom_ptr->ringContactAtom != pdb_j->atom_ptr->atomnumber) &&
					     (pdb_j->atom_ptr->ringContactAtom != pdb_i->atom_ptr->atomnumber))
					    { 
					      /* atoms not bonded to the same atom: >=1_4 interaction */
					      if((pdb_i->atom_ptr->contactAtomNumber != pdb_j->atom_ptr->contactAtomNumber) &&
						 (pdb_i->atom_ptr->twobondAtomNumber != pdb_j->atom_ptr->atomnumber) &&
						 (pdb_j->atom_ptr->twobondAtomNumber != pdb_i->atom_ptr->atomnumber) &&
						 (pdb_i->atom_ptr->ringTwoBondAtom != pdb_j->atom_ptr->atomnumber) &&
						 (pdb_j->atom_ptr->ringTwoBondAtom != pdb_i->atom_ptr->atomnumber))
						{	/*  >=1_4 interaction */ 
						  if((pdb_i->atom_ptr->dihedAtomNumber == pdb_j->atom_ptr->atomnumber) ||
						     (pdb_j->atom_ptr->dihedAtomNumber == pdb_i->atom_ptr->atomnumber) ||
						     (pdb_i->atom_ptr->ringDihedAtom == pdb_j->atom_ptr->atomnumber) ||
						     (pdb_j->atom_ptr->ringDihedAtom == pdb_i->atom_ptr->atomnumber) ||
						     (pdb_i->atom_ptr->contactAtomNumber == pdb_j->atom_ptr->twobondAtomNumber) ||
						     (pdb_j->atom_ptr->contactAtomNumber == pdb_i->atom_ptr->twobondAtomNumber))
						      {	    /* 1_4 interaction */
							 connection_flag=INTRA_ROTAMER_FLAG_14; 
						        /* 1,4 interaction through a ring closure */
						      if((pdb_i->atom_ptr->ringDihedAtom == pdb_j->atom_ptr->atomnumber) || 
							 (pdb_j->atom_ptr->ringDihedAtom == pdb_i->atom_ptr->atomnumber))
							{
							  if(pdb_i->atom_ptr->ringDihedAtom == pdb_j->atom_ptr->atomnumber)
							    {
							      *firstAtom = pdb_i;
							      *dihedAtom = pdb_j;
							    }
							  else
							    {
							      *firstAtom = pdb_j;
							      *dihedAtom = pdb_i;
							    }   
							  
							  if((*firstAtom)->atom_ptr->ringContactAtom!=8)
							  {
							    if((*firstAtom)->atom_ptr->ringContactAtom<(*firstAtom)->atom_ptr->atomnumber)
								step=-1;
							    else
								step=1;
							    *contactAtom = *firstAtom;
							    find_sidechain_atom((*contactAtom),(*firstAtom)->atom_ptr->ringContactAtom, step); 
							  }
							  else
							  {
							      *contactAtom = *firstAtom;
							      stop=0;
							      while((*contactAtom)->seq_position!=0 && stop==0)
							      {
								if((*contactAtom)->atom_ptr->atomnumber==8)
								{
								  if((*contactAtom)->seq_position==(*firstAtom)->seq_position)
								    stop=1;
								  else
								    --*contactAtom; 
								}
								else
								   --*contactAtom;     
							      }
							      if((*contactAtom)->seq_position==0)
							      {
								  *contactAtom = *firstAtom;
								  stop=0;
								  while((*contactAtom)->seq_position!=ENDFLAG && stop==0)
								  {
								    if((*contactAtom)->atom_ptr->atomnumber==8)
								    {
									if((*contactAtom)->seq_position==(*firstAtom)->seq_position)
									    stop=1;
									else
									    ++*contactAtom;     
								    }
								    else
									++*contactAtom;  
								  }	
							      }
							  }
							  
							  if((*firstAtom)->atom_ptr->ringTwoBondAtom!=8)
							  {
							    if((*firstAtom)->atom_ptr->ringTwoBondAtom<(*firstAtom)->atom_ptr->atomnumber)
								step=-1;
							    else
								step=1;
							    *twobondAtom = *firstAtom;
							    find_sidechain_atom((*twobondAtom),(*firstAtom)->atom_ptr->ringTwoBondAtom, step); 
							  }
							  else
							  {
							      *twobondAtom = *firstAtom;
							      stop=0;
							      while((*twobondAtom)->seq_position!=0 && stop==0)
							      {
								if((*twobondAtom)->atom_ptr->atomnumber==8)
								{
								  if((*twobondAtom)->seq_position==(*firstAtom)->seq_position)
								    stop=1;
								  else
								    --*twobondAtom; 
								}
								else
								   --*twobondAtom;     
							      }
							      if((*twobondAtom)->seq_position==0)
							      {
								  *twobondAtom = *firstAtom;
								  stop=0;
								  while((*twobondAtom)->seq_position!=ENDFLAG && stop==0)
								  {
								    if((*twobondAtom)->atom_ptr->atomnumber==8)
								    {
									if((*twobondAtom)->seq_position==(*firstAtom)->seq_position)
									    stop=1;
									else
									   ++*twobondAtom;       
								    }
								    else
									++*twobondAtom;  
								  }	
							      }
							  }
							    
							}
						      else
							{	
							  if((pdb_i->atom_ptr->dihedAtomNumber == pdb_j->atom_ptr->atomnumber) ||
							     (pdb_j->atom_ptr->dihedAtomNumber == pdb_i->atom_ptr->atomnumber))
							    {
							      if(pdb_i->atom_ptr->dihedAtomNumber == pdb_j->atom_ptr->atomnumber)
								{
								  *firstAtom = pdb_i;
								  *dihedAtom = pdb_j;
								}
							      else
								{
								  *firstAtom = pdb_j;
								  *dihedAtom = pdb_i;
								}    
							    }
							  else
							    if((pdb_i->atom_ptr->contactAtomNumber == pdb_j->atom_ptr->twobondAtomNumber) ||
							       (pdb_j->atom_ptr->contactAtomNumber == pdb_i->atom_ptr->twobondAtomNumber))
							      {
								if(pdb_i->atom_ptr->contactAtomNumber == pdb_j->atom_ptr->twobondAtomNumber)
								  {
								    *firstAtom = pdb_j;
								    *dihedAtom = pdb_i;
								  }
								else
								  {
								    *firstAtom = pdb_i;
								    *dihedAtom = pdb_j;
								  }
							      } 
							  
							  if((*firstAtom)->atom_ptr->contactAtomNumber!=8)
							  {    
							    if((*firstAtom)->atom_ptr->contactAtomNumber<(*firstAtom)->atom_ptr->atomnumber)
								step=-1;
							    else
								step=1;
							    *contactAtom = *firstAtom;
							    find_sidechain_atom((*contactAtom),(*firstAtom)->atom_ptr->contactAtomNumber, step); 
							  }
							  else
							  {
							      *contactAtom = *firstAtom;
							      stop=0;
							      while((*contactAtom)->seq_position!=0 && stop==0)
							      {
								if((*contactAtom)->atom_ptr->atomnumber==8)
								{
								  if((*contactAtom)->seq_position==(*firstAtom)->seq_position)
								    stop=1;
								  else
								    --*contactAtom;       
								}
								else
								   --*contactAtom;     
							      }
							      if((*contactAtom)->seq_position==0)
							      {
								  *contactAtom = *firstAtom;
								  stop=0;
								  while((*contactAtom)->seq_position!=ENDFLAG && stop==0)
								  {
								    if((*contactAtom)->atom_ptr->atomnumber==8)
								    {
									if((*contactAtom)->seq_position==(*firstAtom)->seq_position)
									    stop=1;
									else
									  ++*contactAtom;        
								    }
								    else
									++*contactAtom;  
								  }	
							      }
							  }
							  
							  if((*firstAtom)->atom_ptr->twobondAtomNumber!=8)
							  {
							    if((*firstAtom)->atom_ptr->twobondAtomNumber<(*firstAtom)->atom_ptr->atomnumber)
								step=-1;
							    else
								step=1;
							    *twobondAtom = *firstAtom;
							    find_sidechain_atom((*twobondAtom), (*firstAtom)->atom_ptr->twobondAtomNumber, step);
							  }
							  else
							  {
							      *twobondAtom = *firstAtom;
							      stop=0;
							      while((*twobondAtom)->seq_position!=0 && stop==0)
							      {
								if((*twobondAtom)->atom_ptr->atomnumber==8)
								{
								  if((*twobondAtom)->seq_position==(*firstAtom)->seq_position)
								    stop=1;
								  else
								    --*twobondAtom;       
								}
								else
								   --*twobondAtom;     
							      }
							      if((*twobondAtom)->seq_position==0)
							      {
								  *twobondAtom = *firstAtom;
								  stop=0;
								  while((*twobondAtom)->seq_position!=ENDFLAG && stop==0)
								  {
								    if((*twobondAtom)->atom_ptr->atomnumber==8)
								    {
									if((*twobondAtom)->seq_position==(*firstAtom)->seq_position)
									    stop=1;
									else
									  ++*twobondAtom;       
								    }
								    else
									++*twobondAtom;  
								  }	
							      }
							  }  
							} 
						      }
						  else /* >1_4 interaction */
						    {  
							   connection_flag=INTRA_ROTAMER_FLAG; 
						    }  
						}
					    }
					  }
					}			      
				    /* side-bkbn */
				     else  
				       {
				       if(pdb_i->atom_ptr->other_info<0 || pdb_j->atom_ptr->other_info<0)
					 {
					    if(pdb_i->atom_ptr->other_info<0)
					      {
						dummy_bkbn = pdb_i;
						dummy_side = pdb_j;
					      }
					    else
					      {
						dummy_side = pdb_i;
						dummy_bkbn = pdb_j;
					      } 
					    /* 1-4 nonbonded interactions w/ backbone */   
					    /* HB or XG atoms NOT HG  also eliminate (XG/HB)-CA 1-3 interactions  */
					    if((strstr(dummy_side->atom_ptr->atomname,"G")!=NULL && 
					        strstr(dummy_side->atom_ptr->atomname,"HG")==NULL) || 
					       (strstr(dummy_side->atom_ptr->atomname,"HB")!=NULL))
					      {
						if(strcmp(dummy_bkbn->atom_ptr->atomname,"N")==0 || 
						   strcmp(dummy_bkbn->atom_ptr->atomname,"C")==0 ||
						   strstr(dummy_bkbn->atom_ptr->atomname,"HA")!=NULL)
						   {
						    connection_flag=1; 
						    *firstAtom = dummy_side;
						    *dihedAtom = dummy_bkbn;
						    
						    *twobondAtom = dummy_bkbn;
						    find_backbone_atom((*twobondAtom), "CA", 1);
						    
						    *contactAtom =  *twobondAtom;
						    find_backbone_atom((*twobondAtom), "CB", 1);
						   } 
						 else
						   if(strcmp(dummy_bkbn->atom_ptr->atom_ptr->torsiontype,"CT")!=0)  /* not CA or CB */
						    connection_flag=2;  
					      }		
					    else   
					      {  
						/* HG or XD atoms NOT HD */
						if((strstr(dummy_side->atom_ptr->atomname,"D")!=0 && 
						    strstr(dummy_side->atom_ptr->atomname,"HD")==0) ||
						    strstr(dummy_side->atom_ptr->atomname,"HG")!=NULL)     
						  {
						    if(strcmp(dummy_bkbn->atom_ptr->atomname,"CA")==0)
						      {  
						     	connection_flag=1;  
							*firstAtom = dummy_side;
							*dihedAtom = dummy_bkbn;
							
							*twobondAtom = dummy_bkbn;
							find_backbone_atom((*twobondAtom), "CB", 1);
							
							if(dummy_side->atom_ptr->ringDihedAtom!=dummy_bkbn->atom_ptr->atomnumber) 
							  {
							    if(dummy_side->atom_ptr->contactAtomNumber<dummy_side->atom_ptr->atomnumber)
							      step=-1;
							    else
							      step=1;	
							    *contactAtom = dummy_side;
							    find_sidechain_atom((*contactAtom),dummy_side->atom_ptr->contactAtomNumber, step); 
							  }
							
							else  /* 1,4 interaction via ring closure */
							  {
							    if(dummy_side->atom_ptr->ringContactAtom<dummy_side->atom_ptr->atomnumber)
							      step=-1;
							    else
							      step=1;
							    *contactAtom = dummy_side;	
							    find_sidechain_atom((*contactAtom),dummy_side->atom_ptr->ringContactAtom, step); 
							    
							  }
						      }
						    else
						      if(strcmp(dummy_bkbn->atom_ptr->atomname,"CB")!=0)
							connection_flag=2;    
						  }
						else
						  {
						    /* HD or XE atoms NOT HE */   
						    if((strstr(dummy_side->atom_ptr->atomname,"E")!=0 && 
						        strstr(dummy_side->atom_ptr->atomname,"HE")==0) ||
						       (strstr(dummy_side->atom_ptr->atomname,"HD")!=NULL))
						      {
							if(strcmp(dummy_bkbn->atom_ptr->atomname,"CB")==0) 
							  { 
						 	    connection_flag=1;  
							    *firstAtom = dummy_side;
							    *dihedAtom = dummy_bkbn;
							    
							    if(dummy_side->atom_ptr->ringDihedAtom!=dummy_bkbn->atom_ptr->atomnumber) 
							      {
								if(dummy_side->atom_ptr->contactAtomNumber<dummy_side->atom_ptr->atomnumber)
								  step=-1;
								else
								  step=1;	
								*contactAtom = dummy_side;
								find_sidechain_atom((*contactAtom),dummy_side->atom_ptr->contactAtomNumber, step);	 
								
								if(dummy_side->atom_ptr->twobondAtomNumber<dummy_side->atom_ptr->atomnumber)
								  step=-1;
								else
								  step=1;	
								*twobondAtom = dummy_side;	 
								find_sidechain_atom((*twobondAtom), dummy_side->atom_ptr->twobondAtomNumber, step);
							      }  
							    
							    else    /*  1,4 interaction via ring closure   */
							      {
								if(dummy_side->atom_ptr->ringContactAtom<dummy_side->atom_ptr->atomnumber)
								  step=-1;
								else
								  step=1;
								*contactAtom = dummy_side;
								find_sidechain_atom((*contactAtom),dummy_side->atom_ptr->ringContactAtom, step);  	
								
								if(dummy_side->atom_ptr->ringTwoBondAtom<dummy_side->atom_ptr->atomnumber)
								  step=-1;
								else
								  step=1;
								*twobondAtom = dummy_side;
								find_sidechain_atom((*twobondAtom),dummy_side->atom_ptr->ringTwoBondAtom, step);
								
							      }  
							  }
							  else
							    connection_flag=2;	
						      } 
						    else     /*  >1-4 interactions w/ backbone */ 
						    { 
						       connection_flag=2;
						       if(strcmp(dummy_bkbn->atom_ptr->atomname,"CB")==0)
							    connection_flag=INTRA_ROTAMER_FLAG;
						    }      
						  } 
					      }
					 }
				       }
				  }
				}


	return(connection_flag);

}

void build_hbond_burial_virtual_atoms(mini_pdbATOM *probe, mini_pdbATOM *pdb)
{
	int i,p, num_lone_pair,k;	
	mini_pdbATOM *donor_hydrogen, *donor_heavy;
	double bondangle, dihedral, bondlength;
	extern ATOMRESPARAM *BHY_ATOMRESPARAM;
	extern double HBOND_PROBE_BONDLENGTH;

	p=1;
	i=1;
	while(pdb[i].seq_position!=ENDFLAG)
	{

	
		if(pdb[i].atom_ptr->atom_ptr->torsiontype[0] != 'U')
		{
		if(pdb[i].atom_ptr->donorflag[0]!=' ') /* a donor hydrogen or acceptor lone pair */
		{		
			if(pdb[i].atom_ptr->donorflag[0]=='H') /* donor hydrogen; build along same vector */
			{			
				find_donor_heavy(donor_heavy, (&pdb[i]));

				bondangle = pdb[i].atom_ptr->bondangle;
				dihedral = pdb[i].atom_ptr->dihedral;
				bondlength = donor_heavy->atom_ptr->atom_ptr->united_radius + HBOND_PROBE_BONDLENGTH;

				probe[p].seq_position = pdb[i].seq_position;
				probe[p].atom_ptr = BHY_ATOMRESPARAM;

				probe[p].born_radius = i;	/* index of the real atom */
				
		probe[p].coord = add_CARTESIAN(donor_heavy->coord, scalar_multiply_CARTESIAN(bondlength, unitvector(pdb[i].coord, donor_heavy->coord)));

				++p;
			}
			else	/* acceptor */
			{
				num_lone_pair = pdb[i].atom_ptr->donorflag[1] - 48;

				/* acceptor and donor */
				if(pdb[i].atom_ptr->donorflag[0]=='B')
				{
					find_donor_hydrogen(donor_hydrogen, (&pdb[i]));
					
					for(k=1;k<=num_lone_pair;++k)
					{
						/* is donor and acceptor, so build first probe opposite the hydrogen */ 
						if(k==1)
						{
							bondangle = donor_hydrogen->atom_ptr->bondangle;
							dihedral = donor_hydrogen->atom_ptr->dihedral + 180.0;
						}
						if(k==2)
						{
							bondangle = 180.0;
							dihedral = 0.0;
						}
										

					bondlength = pdb[i].atom_ptr->atom_ptr->united_radius + HBOND_PROBE_BONDLENGTH;

					
					probe[p].seq_position = pdb[i].seq_position;
					probe[p].atom_ptr = BHY_ATOMRESPARAM;
					probe[p].coord = pdb[i].coord;

					/*
					contactAtom = &pdb[i];

					twobondAtom = contactAtom;
					find_sidechain_atom(twobondAtom, contactAtom->atom_ptr->contactAtomNumber, -1);
	
					dihedAtom = twobondAtom;	
					find_sidechain_atom(dihedAtom, contactAtom->atom_ptr->twobondAtomNumber, -1);
				
					probe[p].coord = chi2xyz(&contactAtom->coord, &twobondAtom->coord, &dihedAtom->coord, 
                             				&bondlength, &bondangle, &dihedral);
						
					*/

					probe[p].born_radius = i;	/* index of the real atom */

					++p;

					}
				}
				else if(pdb[i].atom_ptr->donorflag[0]=='A') /* just an acceptor */
				{
					for(k=1;k<=num_lone_pair;++k)
					{
						if(k==1)
						{	
							bondangle = 109.5;
							dihedral = 0.0;
						}
						if(k==2)
						{	
							bondangle = 109.5;
							dihedral = 180.0;
						}
						

						bondlength = pdb[i].atom_ptr->atom_ptr->united_radius + HBOND_PROBE_BONDLENGTH;

						probe[p].seq_position = pdb[i].seq_position;
						probe[p].atom_ptr = BHY_ATOMRESPARAM;
						probe[p].coord = pdb[i].coord;

						/*
						contactAtom = &pdb[i];

						twobondAtom = contactAtom;
						find_sidechain_atom(twobondAtom, contactAtom->atom_ptr->contactAtomNumber, -1);
		
						dihedAtom = twobondAtom;	
						find_sidechain_atom(dihedAtom, contactAtom->atom_ptr->twobondAtomNumber, -1);
				
						probe[p].coord = chi2xyz(&contactAtom->coord, &twobondAtom->coord, &dihedAtom->coord, 
                             				&bondlength, &bondangle, &dihedral);
						*/

						probe[p].born_radius = i;	/* index of the real atom */

					++p;

					}
				}
			}
		}
		}
		++i;
	}
	probe[p] = pdb[i];

}


/* returns sasa of solvent-accessible hbond probe atoms */
double exposed_hbond_groups(mini_pdbATOM *pdb)
{
  double exposed_hbond;

  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, *hbond_donor;
  int		box_atoms;		/* Total atoms in area calculation */
  static mini_pdbATOM *hbond_probe_pdb=NULL;
  double sasa_pts=0;
  static SASA_ATOM *atm=NULL;		/* Atom array */
  static int max_atoms=0;
  extern SASA_ATOM   **box;
  extern int BHY_CLASS;
  extern int *BOX_OFFSET;
  extern unsigned char  SASA_CLOSEST_POINT[][200];
  extern SASA_point *point;
  extern double **SASA_RADII2_DIFF, ATOM_POINTS;
  extern int  sasa_count[], atm_pts_minus_MAXPTS;
	  

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


  pdb0 = pdb;
  ++pdb;
  box_atoms=0;
  while(pdb->seq_position!=ENDFLAG)
    {
      /* 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;

  pdb0 = hbond_probe_pdb;
  ++hbond_probe_pdb;
  while(hbond_probe_pdb->seq_position!=ENDFLAG)
    {
      hbond_probe_pdb->sasa=0;
      
	  atm[box_atoms].pdb = hbond_probe_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;	
	
      ++hbond_probe_pdb;
    }
    hbond_probe_pdb = pdb0;

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

  /* 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];
		      
		      /* one and only one of the atoms must be of type BHY */
			 if(a1->pdb->atom_ptr->atom_ptr->atom_class == BHY_CLASS ||
			 a2->pdb->atom_ptr->atom_ptr->atom_class == BHY_CLASS &&
			 a1->pdb->atom_ptr->atom_ptr->atom_class != a2->pdb->atom_ptr->atom_ptr->atom_class)
		      {
		      
		      /* 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 only */
				  if(a1->pdb->atom_ptr->atom_ptr->atom_class == BHY_CLASS) /* a1 is a BHY */
				    {
				      /* 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++));
				    }
				  else	/* a2 is a BHY */
				    {
				      /* 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->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 BHY atoms */
  for (i=0;i<box_atoms;i++)
    {
      /* Count total number of exposed surface points */

	if(atm[i].pdb->atom_ptr->atom_ptr->atom_class == BHY_CLASS)
	{
      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;

	}
    }		



   j=1;	  
   while(hbond_probe_pdb[j].seq_position != ENDFLAG)
   {
	i = (int)hbond_probe_pdb[j].born_radius;

	if(pdb[i].atom_ptr->donorflag[0]=='H')	// hydrogen
	{
		find_donor_heavy(hbond_donor, (&pdb[i]));  // donor

		if(hbond_probe_pdb[j].sasa > 0)	
		{
				pdb[i].hbond_satisfied=1;
				++exposed_hbond;

				
				hbond_donor->hbond_satisfied=1;
				++exposed_hbond;
		}
		else if(hbond_donor->sasa>0)
		{
			hbond_donor->hbond_satisfied=1;
			++exposed_hbond;
		}
		
	}	
	else 	// acceptor
	{
		
		if(pdb[i].sasa > 0)	
		{
			++exposed_hbond;
			pdb[i].hbond_satisfied=1;
		}
		
	}	
	 

	++j;
   }

  return(exposed_hbond);
}


ENERGY energy_for_atom_of_interest(mini_pdbATOM *pdb, int index_of_atom_of_interest, int electr_flag, int sasa_flag, int torsion_flag) 
{
  double r2; /* squared distance between two atoms */
  double EE, r;
  ENERGY energy;
  double E_coulomb, E_coulomb_14, E_vdw_14, E_torsion;
  double  alpha, f;
  double xmax, ymax, zmax, xmin, ymin, zmin;
  int flag2, j_ctr;
  mini_pdbATOM *firstAtom, *dihedAtom, *twobondAtom, *contactAtom;
  COULOMBIC *first_coulombic, *coulombic;
  double dx, dy, dz;
  VDW_ENERGY vdwE;
  extern double **CHARGE_PRODUCT;
  mini_pdbATOM *atom_of_interest;

    
  ENERGY_0(energy); VDW_ENERGY_0(vdwE);
  E_coulomb=0;  E_coulomb_14=0; E_vdw_14=0; E_torsion=0;  
  firstAtom = NULL; contactAtom = NULL; twobondAtom = NULL; dihedAtom = NULL;
  first_coulombic = NULL; coulombic = NULL; 

  atom_of_interest = &pdb[index_of_atom_of_interest];

   if(atom_of_interest->seq_position==ENDFLAG) 
    	return(energy);
    if(atom_of_interest->atom_ptr->atom_ptr->atom_class==ENDFLAG)  
    	return(energy);

	/* initialize COULOMBIC linked-list */

    if(electr_flag!=0)
    {
    coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC));
    first_coulombic = coulombic;   
    coulombic->next = NULL;  
    }
  
	  xmax=atom_of_interest->coord.x+FORCEFIELD_DISTANCE_CUTOFF;
	  xmin=atom_of_interest->coord.x-FORCEFIELD_DISTANCE_CUTOFF;
	  ymax=atom_of_interest->coord.y+FORCEFIELD_DISTANCE_CUTOFF;
	  ymin=atom_of_interest->coord.y-FORCEFIELD_DISTANCE_CUTOFF;
	  zmax=atom_of_interest->coord.z+FORCEFIELD_DISTANCE_CUTOFF;
	  zmin=atom_of_interest->coord.z-FORCEFIELD_DISTANCE_CUTOFF;
	  
      j_ctr= 1; 
      while(pdb[j_ctr].seq_position!=ENDFLAG)  /* j loop */
	{
	if(pdb[j_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG)
	{
	  /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF
	     centered at pdb[j_ctr].coord */
	  if(pdb[j_ctr].coord.x<=xmax)
	    if(pdb[j_ctr].coord.x>=xmin)
	      if(pdb[j_ctr].coord.y<=ymax)
		if(pdb[j_ctr].coord.y>=ymin)
		  if(pdb[j_ctr].coord.z<=zmax)
		    if(pdb[j_ctr].coord.z>=zmin)
		      {
			dx = pdb[j_ctr].coord.x - atom_of_interest->coord.x;
			dy = pdb[j_ctr].coord.y - atom_of_interest->coord.y;
			dz = pdb[j_ctr].coord.z - atom_of_interest->coord.z;
			r2 = dx*dx + dy*dy + dz*dz;  
			if(r2>=TINY &&  r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD)  /* FORCEFIELD_DISTANCE_CUTOFF cutoff */  
			  {
			    flag2=0;  
			    r=sqrt(r2);
			      
					/* get the connectivity between the atoms; if 1,4, also get pointers to the dihedral atoms */
				flag2 = get_atom_connectivity(r2, &pdb[index_of_atom_of_interest], &pdb[j_ctr], 
												&firstAtom, &contactAtom, &twobondAtom, &dihedAtom);

			/* calculate desolvation for born radii calcs  */
			     if(atom_of_interest->atom_ptr->charge!=0) /*  i charged  */
			      {   
			        
				if(pdb[j_ctr].atom_ptr->charge!=0)  /* i and j charged */
				  {
				     
				    if(flag2<=3)	/* save distances for born calcs later */
				       {
					    
						if(electr_flag!=0)
						{
						coulombic->index_i =  index_of_atom_of_interest;
						coulombic->index_j = j_ctr; 
						coulombic->r2 = r2;
						coulombic->flag = flag2;
					    
						coulombic->next = (COULOMBIC *)malloc(sizeof(COULOMBIC));  /* allocate memory for next element */
						coulombic = coulombic->next;  
						coulombic->next = NULL; /* last element unless replaced */
					    }

					  } 
				    }  		     
			      }  	
			    
			     if(flag2!=0) 
			    {
				  if(flag2==2)    /*   >1,4 interaction  */
				    {
						E_coulomb += coulomb_energy_q2_per_angstrom(r, atom_of_interest->atom_ptr, pdb[j_ctr].atom_ptr);
						add_atompair_vdw_energy(&vdwE, atom_of_interest->atom_ptr, pdb[j_ctr].atom_ptr, r2, r);

						if(electr_flag == 1)
							energy.E_hbond += hbond_energy(&pdb[index_of_atom_of_interest], &pdb[j_ctr], r, r2);

				    } 
				   if(flag2==1)       /* 1,4 */
					if(torsion_flag==1)
				    {
					
					  E_vdw_14 += vdw_energy_14(atom_of_interest->atom_ptr, pdb[j_ctr].atom_ptr, r2);
				      E_coulomb_14 += coulomb_energy_q2_per_angstrom(r, atom_of_interest->atom_ptr, pdb[j_ctr].atom_ptr);
				      E_torsion += torsion_energy(firstAtom, contactAtom, twobondAtom, dihedAtom); 
				     }   
			    }   
			}
		   }
	}    

	++j_ctr;
	}
  
    /* calculate born energies */
  energy.E_born = 0; energy.E_pol=0;
  if(atom_of_interest->atom_ptr->charge!=0)
  if(electr_flag!=0)
    {
      
	/* born cross-polarization energies; use the coulombic linked-list with the born radii calculated above */
      coulombic = first_coulombic;
      while(coulombic->next!=NULL)
	{
	  alpha = atom_of_interest->born_radius*pdb[coulombic->index_j].born_radius;
	   
	  f = f_of_r2_alpha(coulombic->r2, alpha);
	  
	  EE = CHARGE_PRODUCT[atom_of_interest->atom_ptr->charge_class][pdb[coulombic->index_j].atom_ptr->charge_class]/f; 
	  
	/* prevent favorable born pair energies from overwhelming lennard-jones for close, nonbonded atoms */
	  if(coulombic->flag == 2)
	  {
	      if(EE > POL_ATTRACTION_CAP)	 
		EE = POL_ATTRACTION_CAP;
	  }

	  energy.E_pol += BORN_CONST(f)*EE;
	
	  coulombic = coulombic->next;      
	}
    
    coulombic = first_coulombic;
    freeCOULOMBIC(coulombic);

    coulombic = NULL; 
    first_coulombic = NULL;  

	energy.E_pol = COULOMB_CONST*energy.E_pol;
	energy.E_born = HALF_COULOMB_CONST*BORN_CONST(atom_of_interest->born_radius)*atom_of_interest->atom_ptr->charge2/atom_of_interest->born_radius;
  }

	/* add up energies, multiply with constants to get the correct units, etc */  

  energy.E_coulomb = COULOMB_CONST_OVER_INTERNAL_DIELECTRIC*E_coulomb;

  energy.E_1_4 = WEIGHT_1_4*(NONBOND_FACTOR_1_4*(E_vdw_14 + COULOMB_CONST*E_coulomb_14) + E_torsion);

  energy.E_vdw = sum_VDW_ENERGY(vdwE);	 

	energy.E_hbond = HBOND_FUNCTION_FLAG*energy.E_hbond;

  if(sasa_flag!=0)
  if(atom_of_interest->atom_ptr->atom_ptr->torsiontype[0]=='C' || 
	atom_of_interest->atom_ptr->atom_ptr->torsiontype[0]=='H' || 
	atom_of_interest->atom_ptr->atom_ptr->torsiontype[0]=='S' &&
	atom_of_interest->atom_ptr->atom_ptr->torsiontype[0]!='U')
  		energy.E_sasa = HYDROPHOBIC_ASP*atom_of_interest->sasa;

   /* divide pair energies by 2 to account for overcounting */
	energy.E_pol = 0.5*energy.E_pol;
	energy.E_coulomb = 0.5*energy.E_coulomb;
	energy.E_1_4 = 0.5*energy.E_1_4;
	energy.E_vdw = 0.5*energy.E_vdw;

	/*
	printf("%d\t%s\t%lf\t%lf\t%lf\n",atom_of_interest->seq_position, atom_of_interest->atom_ptr->atomname,
			energy.E_vdw , 	energy.E_1_4, energy.E_coulomb + energy.E_pol + energy.E_born);
	*/

 return energy;
  
}