/* EGAD: ligand_stuff.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 reading ligand parameter files and for generating and manipulating 
	ligand structures.

*/


#include "ligand_stuff.h"

double *STATIONARY_LIGAND_TRANSFORM_MATRIX;

/* Extracts the transform matrix for a ligand given its coords in pdb. 
	resparam must be assigned to the proper ligand by the calling function.
	The matrix is read into the chi array (allocated by calling function).
	If the ligand "rotamers" are input directly as coordinates, this function places the 
		coords in pdb into the last "rotamer" for this ligand.
	This function is called by extract_dihedral; rotatable bond dihedrals are calculated there if needed. */
void extract_transform_matrix_for_ligand(pdbATOM *pdb, double *chi, RESPARAM *resparam)
{
	int i,j,k,chi_index,endSideAtomNumber, previous_atoms[6];
	int a,b,c,stopflag;
	double **matrix_A, *new_dimensions, *soln_vector;
	CARTESIAN *new_coordinates, *base_coordinates;
	extern int MAX_RES_SIZE;

	if(resparam->ligand_flag==0)
	{
		fprintf(stderr,"ERROR extract_transform_matrix_for_ligand is being called with non-ligand %s\n",resparam->residuetype);
		fprintf(stderr,"dumping core.....\n");
		abort();
	}

	
	/* 'rotamers' are coordinates */
	
	if(resparam->rotamerlib_ptr->numChi==1) 
	{
		j=0; i=9; endSideAtomNumber = resparam->numberofAtoms+1;
	
			while(strcmp(resparam->atom[i].atomtype,"U")==0)	/* move past virtual atoms if necessary */
				++i;


		chi[1] = resparam->rotamerlib_ptr->number_of_rotamers;
		
		if(resparam->rotamerlib_ptr->rotamer[resparam->rotamerlib_ptr->number_of_rotamers].ligand_coords == NULL)
			resparam->rotamerlib_ptr->rotamer[resparam->rotamerlib_ptr->number_of_rotamers].ligand_coords = 
				(pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM));

			while(i<=endSideAtomNumber)
			{
			++j;
			k=1;
			while(strcmp(pdb[k].atomname,resparam->atom[i].atomname)!=0 && pdb[k].seq_position!=ENDFLAG)
				++k;
			if(pdb[k].seq_position==ENDFLAG)
			{
				fprintf(stderr,"ERROR Cannot find %s in the coords or ligamer for ligand %s\n",
						resparam->atom[i].atomname,resparam->residuetype);
				exit(1);
			}
			resparam->rotamerlib_ptr->rotamer[resparam->rotamerlib_ptr->number_of_rotamers].ligand_coords[j] = pdb[k];
		
			++i;
			}
			++j;
			resparam->rotamerlib_ptr->rotamer[resparam->rotamerlib_ptr->number_of_rotamers].ligand_coords[j].seq_position=ENDFLAG;
			resparam->rotamerlib_ptr->rotamer[resparam->rotamerlib_ptr->number_of_rotamers].ligand_coords[j].atom_ptr=NULL;			
			strcpy(resparam->rotamerlib_ptr->rotamer[resparam->rotamerlib_ptr->number_of_rotamers].ligand_coords[j].residuetype,"END");

		return;
	}

	/* new_dimension[i] = a*x_i + b*y_i + c*z_i + d */
	/* given new coordinates and base coordinates, one needs only 4 points i to solve 
		a 4 eqn. system for the 4 transformation matrix elements a->d for this dimension */
	/* this assumes fixed internal geometry between these four points (preceeding the 1st rotatable bond) */
	/* rounding error can lead to scaling and other errors!!!!!!!! */


	if(resparam->rotamerlib_ptr->ligand_base_coords==NULL)
	{
		fprintf(stderr,"ERROR ligand_base_coords must be defined in extract_transform_matrix_for_ligand for ligand %s\n",
				resparam->residuetype);
		fprintf(stderr,"dumping core.....\n");
		abort();
	}

	new_coordinates = (CARTESIAN *)calloc(5,sizeof(CARTESIAN));
	base_coordinates = (CARTESIAN *)calloc(5,sizeof(CARTESIAN));

	new_dimensions = (double *)calloc(5,sizeof(double));
	soln_vector = (double *)calloc(5,sizeof(double));
	matrix_A = (double **)calloc(5,sizeof(double *));
	for(i=1;i<=4;++i)
	{
		matrix_A[i]	= (double *)calloc(5,sizeof(double));
		previous_atoms[i]=ENDFLAG;
	}
	previous_atoms[5]=ENDFLAG;

	j=0; i=9; endSideAtomNumber = resparam->numberofAtoms+1;
	
	while(strcmp(resparam->atom[i].atomtype,"U")==0)	/* move past virtual atoms if necessary */
		++i;

	while(j<4)
	{
		if(resparam->atom[i].dihedral>=0)
		{
			++j;
			k=1;
			while(strcmp(pdb[k].atomname,resparam->atom[i].atomname)!=0 && pdb[k].seq_position!=ENDFLAG)
				++k;
			if(pdb[k].seq_position==ENDFLAG)
			{
				fprintf(stderr,"ERROR Cannot find %s in the coords or ligamer for ligand %s\n",
						resparam->atom[i].atomname,resparam->residuetype);
				exit(1);
			}
			new_coordinates[j] = pdb[k].coord;
			k=1;
			while(strcmp(resparam->rotamerlib_ptr->ligand_base_coords[k].atomname,resparam->atom[i].atomname)!=0)
				++k;
			base_coordinates[j] = resparam->rotamerlib_ptr->ligand_base_coords[k].coord;
			previous_atoms[j]=resparam->atom[i].atomnumber;
			++i;
		}
		else
			{
				
				// variable dihedral for this atom....find another atom that is dependent only on the preceeding rigid atoms
				stopflag=0;
				while(stopflag==0)
				{
					++i;	
					
					if(strcmp(resparam->atom[i].atomname,"zzz")==0)
					{
						fprintf(stderr,"ERROR Need at least 4 atoms preceeding the first rotatable bond for %s\n",resparam->residuetype);
						exit(1);
					}

					if(resparam->atom[i].dihedral>=0)	// fixed-dihedral for this atom...is it dependent only on the previous atoms?
					{
						k=1; a=0; b=0; c=0;
						while(previous_atoms[k]!=ENDFLAG)
						{
							if(resparam->atom[i].contactAtomNumber == previous_atoms[k] || 
								resparam->atom[i].contactAtomNumber < previous_atoms[1])
								a=1;
							if(resparam->atom[i].twobondAtomNumber == previous_atoms[k] || 
								resparam->atom[i].twobondAtomNumber < previous_atoms[1])
								b=1;
							if(resparam->atom[i].dihedAtomNumber == previous_atoms[k] || 
								resparam->atom[i].dihedAtomNumber < previous_atoms[1])
								c=1;
							++k;
						}
						if(a==1 && b==1 && c==1)
							stopflag=1;
					}
				}
			}
		
	}	
	
	for(i=1;i<=4;++i)
	{
		matrix_A[i][1] = base_coordinates[i].x;
		matrix_A[i][2] = base_coordinates[i].y;
		matrix_A[i][3] = base_coordinates[i].z;
		matrix_A[i][4] = 1;
	}

	chi_index=1;
	for(k=1;k<=3;++k)
	{
		for(i=1;i<=4;++i)
		{
			if(k==1)
				new_dimensions[i] = new_coordinates[i].x;
			else if(k==2)
				new_dimensions[i] = new_coordinates[i].y;
			else if(k==3)
				new_dimensions[i] = new_coordinates[i].z;
		}

		cramers_rule(4, matrix_A, soln_vector, new_dimensions);

		for(i=1;i<=4;++i)
		{
			chi[chi_index] = soln_vector[i]; 
			++chi_index;
		}

	}

	free_memory(soln_vector); free_memory(new_dimensions); 
	for(i=1;i<=4;++i)
		free_memory(matrix_A[i]);
	free_memory(matrix_A);
	free_memory(new_coordinates); free_memory(base_coordinates);
	
}


/* called by build_a_Sidechain.
	This function loads the SIDECHAIN ligand with coords corresponding to "rotamer" chi.
	For "rotamers" = coordinates, chi = rotamerlib index.
	For "rotamers" = transform matrix + dihedrals, chi = matrix + dihedrals.
	resparam must be assigned to the proper ligand by the calling function */

void build_a_Ligand(SIDECHAIN *ligand, RESPARAM *resparam, double chi[])
{
	int i,n,k,a,b,c,sideAtomNumber, endSideAtomNumber, startSideAtomNumber, *flag, *rigid_atoms, movable_flag;	
	double *real_chi;
	double dihedral, dihed2;
	CARTESIAN *dum;
	

	if(resparam->ligand_flag==0)
	{
		fprintf(stderr,"ERROR build_a_Ligand is being called with non-ligand %s\n",resparam->residuetype);
		fprintf(stderr,"dumping core.....\n");
		abort();
	}


	if(resparam->rotamerlib_ptr->numChi==1) /* use inputted coordinates directly */
	{	
			i=1;sideAtomNumber=9;

			while(strcmp(resparam->atom[sideAtomNumber].atomtype,"U")==0)
			{
				CARTESIAN_zero(ligand->atom[sideAtomNumber].coord);
				ligand->atom[sideAtomNumber].atom_ptr = &resparam->atom[sideAtomNumber];
				++sideAtomNumber;
			}

			while(resparam->rotamerlib_ptr->rotamer[(int)chi[1]].ligand_coords[i].seq_position!=ENDFLAG)
			{
				ligand->atom[sideAtomNumber].coord = 
					resparam->rotamerlib_ptr->rotamer[(int)chi[1]].ligand_coords[i].coord;
				ligand->atom[sideAtomNumber].atom_ptr = &resparam->atom[sideAtomNumber];

				++sideAtomNumber; ++i;
			}
			ligand->atom[sideAtomNumber].atom_ptr = NULL;

			return;
	}

		/* "rotamers" are transform matricies; transform from base coords  */

			i=1;sideAtomNumber=9;

			while(strcmp(resparam->atom[sideAtomNumber].atomtype,"U")==0)
			{
				CARTESIAN_zero(ligand->atom[sideAtomNumber].coord);
				ligand->atom[sideAtomNumber].atom_ptr = &resparam->atom[sideAtomNumber];
				++sideAtomNumber;
			}

			while(resparam->rotamerlib_ptr->ligand_base_coords[i].seq_position!=ENDFLAG)
			{
				ligand->atom[sideAtomNumber].coord = 
					tranform_coordinates(resparam->rotamerlib_ptr->ligand_base_coords[i].coord, chi);
				ligand->atom[sideAtomNumber].atom_ptr = &resparam->atom[sideAtomNumber];

				++sideAtomNumber; ++i;
			}
			ligand->atom[sideAtomNumber].atom_ptr = NULL;

		if(resparam->rotamerlib_ptr->numChi > 12)  /* rotatable bonds too */
		{
			real_chi = (double *)calloc(resparam->num_rotatable_bonds+2,sizeof(double));
			flag = (int *)calloc(MAX_RES_SIZE,sizeof(int));
			rigid_atoms = (int *)calloc(MAX_RES_SIZE,sizeof(int));
			dum = (CARTESIAN *)calloc(MAX_RES_SIZE,sizeof(CARTESIAN));
			
			k=1;
			while(k<MAX_RES_SIZE)
			{
				rigid_atoms[k]=ENDFLAG;
				++k;
			}

			n=1;
			for(i=13;i<=resparam->rotamerlib_ptr->numChi;++i)
			{
				real_chi[n] = chi[i];	
				++n;
			}
			real_chi[n]=0;

			
			endSideAtomNumber = resparam->numberofAtoms+1;

			sideAtomNumber=9;

			while(strcmp(resparam->atom[sideAtomNumber].atomtype,"U")==0)
			{
				CARTESIAN_zero(ligand->atom[sideAtomNumber].coord);
				ligand->atom[sideAtomNumber].atom_ptr = &resparam->atom[sideAtomNumber];
				++sideAtomNumber;
			}

			k=1;
			while(resparam->atom[sideAtomNumber].dihedral>=0)
			{
				dum[sideAtomNumber] = ligand->atom[sideAtomNumber].coord;
				rigid_atoms[k] = sideAtomNumber;

				/*
				printf("start\t%s\t\t%lf\t%lf\t%lf\n",resparam->atom[sideAtomNumber].atomname,
					dum[sideAtomNumber].x,dum[sideAtomNumber].y,dum[sideAtomNumber].z);
				*/

				++k;
				++sideAtomNumber;
			}
			startSideAtomNumber = sideAtomNumber; /* index of atom w/ first rotatable bond */
			
			while(strcmp(resparam->atom[sideAtomNumber].atomtype,"zzz")!=0)
			{
				dum[sideAtomNumber] = ligand->atom[sideAtomNumber].coord;
				++sideAtomNumber;
			}	

			n=0;
  			for(sideAtomNumber=startSideAtomNumber;sideAtomNumber<=endSideAtomNumber;++sideAtomNumber)
    			{

			/*
				printf("here\t%s\t%lf\t\t%lf\t%lf\t%lf\n",resparam->atom[sideAtomNumber].atomname,resparam->atom[sideAtomNumber].dihedral,dum[sideAtomNumber].x,dum[sideAtomNumber].y,dum[sideAtomNumber].z);
			*/

      				if(resparam->atom[sideAtomNumber].dihedral<0)   /* variable chi angle */
				{
	  				++n; flag[sideAtomNumber]=0;
	  				dihedral=real_chi[n]; 
	  				/* change dihedrals of other atoms that will change with this one */
	  				for(i=startSideAtomNumber;i<=endSideAtomNumber;++i) 
	    					if(resparam->atom[i].contactAtomNumber==resparam->atom[sideAtomNumber].contactAtomNumber && 
	       						i!=sideAtomNumber)
	      					{
							flag[i]=1; /* mark those atoms w/ changed dihedrals so they won't be changed when sideAtomNumber=i */
 							dihed2=dihedral+(resparam->atom[i].dihedral -(-1)*(resparam->atom[sideAtomNumber].dihedral));
							dum[i] = chi2xyz(&dum[resparam->atom[i].contactAtomNumber],
											&dum[resparam->atom[i].twobondAtomNumber],
											&dum[resparam->atom[i].dihedAtomNumber],
											&(resparam->atom[i].bondlength),
											&(resparam->atom[i].bondangle),
											&dihed2);


	      					}
				}
      				else
	  				dihedral = resparam->atom[sideAtomNumber].dihedral;
      				if(flag[sideAtomNumber]!=1) /* won't calc. coords for those atoms we already took care of above */ 
				{
					// don't calc new coords for rigid atoms ie: dependent only on rigid atoms
					movable_flag=1; 	// assume moving, unless actually rigid
					
					if(resparam->atom[sideAtomNumber].dihedral>0)
					{
					k=1; a=0; b=0; c=0;
					while(rigid_atoms[k]!=ENDFLAG)
					{					
						if(resparam->atom[sideAtomNumber].contactAtomNumber == rigid_atoms[k] ||
							resparam->atom[sideAtomNumber].contactAtomNumber < rigid_atoms[1])
							a=1;
						if(resparam->atom[sideAtomNumber].twobondAtomNumber == rigid_atoms[k] ||
							resparam->atom[sideAtomNumber].twobondAtomNumber < rigid_atoms[1])
							b=1;
						if(resparam->atom[sideAtomNumber].dihedAtomNumber == rigid_atoms[k] ||
							resparam->atom[sideAtomNumber].dihedAtomNumber < rigid_atoms[1])
							c=1;
						
						if(sideAtomNumber == rigid_atoms[k])
						{
							a=1; b=1; c=1;
						}

						++k;
					}
					if(a==1 && b==1 && c==1)	// this is a rigid atom
					{
						movable_flag=0;
						rigid_atoms[k]=sideAtomNumber;

					}
					}

					if(movable_flag==1)
						dum[sideAtomNumber] = chi2xyz(&dum[resparam->atom[sideAtomNumber].contactAtomNumber],
				      								&dum[resparam->atom[sideAtomNumber].twobondAtomNumber],
				      								&dum[resparam->atom[sideAtomNumber].dihedAtomNumber],
				      								&(resparam->atom[sideAtomNumber].bondlength),
				      								&(resparam->atom[sideAtomNumber].bondangle),
				      								&dihedral);
				}
    			}

			ligand->atom[endSideAtomNumber+1].atom_ptr = NULL;
			for(sideAtomNumber=startSideAtomNumber;sideAtomNumber<=endSideAtomNumber;++sideAtomNumber)
			{
				ligand->atom[sideAtomNumber].coord = dum[sideAtomNumber];
			
				/*				
				printf("%s\t%lf\t%lf\t%lf\n",resparam->atom[sideAtomNumber].atomname,ligand->atom[sideAtomNumber].coord.x,ligand->atom[sideAtomNumber].coord.y,ligand->atom[sideAtomNumber].coord.z);
				*/
			}

			free_memory(real_chi); free_memory(flag); free_memory(dum); free_memory(rigid_atoms);

		}

}


/* called by read_forcefield.cpp: read_rotamers if GENERATE_LIGAMERS is defined */
void generate_ligamers(FILE *input_file_ptr, ROTAMERLIB *rotamerlib, char *line)
{
	char *dummyfilename;
	FILE *file_ptr;
	double dx, dy, dz, theta;
	double *translate_rotate_array, *transform_matrix;
	double **rotation_axis;
	int a,i;	
	double rotation_stepsize, translation_stepsize, max_displacement;
	extern int GET_PID;
	
	rotation_stepsize = 20.0; translation_stepsize = 1.0;  max_displacement = 2.0;

	dummyfilename = (char *)calloc(MAXLINE,sizeof(char));

	if(word_count(line)==5)
		sscanf(line,"%s %s %lf %lf %lf",rotamerlib->residuetype, dummyfilename, 
						&rotation_stepsize, &translation_stepsize, &max_displacement);
	
	rotamerlib->generate_ligamers_flag=1;

	transform_matrix = (double *)calloc(14,sizeof(double));
	translate_rotate_array = (double *)calloc(9,sizeof(double));

	translate_rotate_array[1] = 0; translate_rotate_array[2] = 0; translate_rotate_array[3] = 0;
	for(i=1;i<=12;++i)
		transform_matrix[i] = 0;

	rotation_axis = (double **)calloc(15,sizeof(double *));
	for(i=1;i<=13;++i)
		rotation_axis[i] = (double *)calloc(5,sizeof(double));

	/* 13 arbitrary axes of rotations; x,y,z axes and all diagonals */
	rotation_axis[1][1] = 1; rotation_axis[1][2] = 1; rotation_axis[1][3] = 1; 

	rotation_axis[2][1] = 1; 	rotation_axis[2][2] = 1; 	rotation_axis[2][3] = 0; 
	rotation_axis[3][1] = 1; 	rotation_axis[3][2] = 0; 	rotation_axis[3][3] = 1; 
	rotation_axis[4][1] = 0; 	rotation_axis[4][2] = 1; 	rotation_axis[4][3] = 1; 

	rotation_axis[5][1] = 1; 	rotation_axis[5][2] = -1; 	rotation_axis[5][3] = 0; 
	rotation_axis[6][1] = 1; 	rotation_axis[6][2] = 0; 	rotation_axis[6][3] = -1; 
	rotation_axis[7][1] = 0; 	rotation_axis[7][2] = 1; 	rotation_axis[7][3] = -1; 

	rotation_axis[8][1] = 1; 	rotation_axis[8][2] = 0; 	rotation_axis[8][3] = 0; 
	rotation_axis[9][1] = 0; 	rotation_axis[9][2] = 1; 	rotation_axis[9][3] = 0; 	
	rotation_axis[10][1] = 0; 	rotation_axis[10][2] = 0; 	rotation_axis[10][3] = 1; 

	rotation_axis[11][1] = 1; 	rotation_axis[11][2] = -1; 	rotation_axis[11][3] = -1; 
	rotation_axis[12][1] = -1; 	rotation_axis[12][2] = 1; 	rotation_axis[12][3] = -1; 
	rotation_axis[13][1] = -1; 	rotation_axis[13][2] = -1; 	rotation_axis[13][3] = 1; 


	sprintf(dummyfilename,"ligamers.%s.%d", rotamerlib->residuetype, GET_PID);
	file_ptr = fopen_file(dummyfilename,"w");

	
	fprintf(file_ptr,"1 0 0 0 0 1 0 0 0 0 1 0\n");  // unchanged base coordinates 

	for(dx = -max_displacement; dx <= max_displacement ; dx += translation_stepsize)
	{

	translate_rotate_array[1] = dx;

	for(dy = -max_displacement; dy <= max_displacement ; dy += translation_stepsize)
	{

	translate_rotate_array[2] = dy;

	for(dz = -max_displacement; dz <= max_displacement ; dz += translation_stepsize)
	{

	translate_rotate_array[3] = dz;

	for(a=1;a<=13;++a)
	{

	translate_rotate_array[4] = rotation_axis[a][1];	
	translate_rotate_array[5] = rotation_axis[a][2];
	translate_rotate_array[6] = rotation_axis[a][3];
	
	for(theta = -180 ; theta < 180 ; theta += rotation_stepsize)
	{
		translate_rotate_array[7] = theta;
	
		translate_rotate_array_to_transform_matrix(translate_rotate_array, transform_matrix);

		check_transform_matrix(transform_matrix);

			for(i=1;i<=12;++i)
				fprintf(file_ptr,"%.14lf ",transform_matrix[i]);
			fprintf(file_ptr,"\n");

	} // end theta


	}  // axis of rotation list

	
	} // end dz
	} // end dy
	} // end dx

	fprintf(file_ptr,"STOP\n");
	fclose(file_ptr);
	
	fclose(input_file_ptr);

	input_file_ptr = fopen_file(dummyfilename,"r");

	for(i=1;i<=13;++i)
		free_memory(rotation_axis[i]);
	free_memory(rotation_axis);

	free_memory(dummyfilename);
	free_memory(translate_rotate_array);
	free_memory(transform_matrix);
	
}


#include "output_stuff.h"

/* read ligand coordinates into ligand_coords (allocated by calling function).
	file_ptr at "START" line of the file, immediately preeceding the first coordinate 
	Reads until "END". */

void read_ligand_coords(FILE *file_ptr, pdbATOM *ligand_coords)
{
	char *line, *dummystring;
	int i,a,b; 
        static ATOMRESPARAM *dummy_atom_ptr=NULL;

	if(dummy_atom_ptr==NULL)
	{
		dummy_atom_ptr=(ATOMRESPARAM *)malloc(sizeof(ATOMRESPARAM));
		dummy_atom_ptr->atom_ptr =(ATOMPARAM *)malloc(sizeof(ATOMPARAM));
		dummy_atom_ptr->atom_ptr->atom_class = ENDFLAG2;
	}

	line = (char *)calloc(MXLINE_INPUT,sizeof(char));
	dummystring = (char *)calloc(25,sizeof(char));

	i=1;
	while(fgets(line,MXLINE_INPUT,file_ptr)!=NULL && strncmp(line,"END",3)!=0)
	{	
	   if(strncmp(line,"ATOM",4)==0 || strncmp(line,"LIG",3)==0)
	   {
		for(a=0;a<20;++a)
                dummystring[a] = '\0';
             a=0; b = 6;
             while(b<=10)
             {   
                dummystring[a] =  line[b];
                ++a; ++b; 
             }
             sscanf(dummystring, "%d", &ligand_coords[i].atom_number);

		 for(a=0;a<20;++a)
                dummystring[a] = '\0';
             a=0; b = 12;
             while(b<=15)
             {   
                dummystring[a] =  line[b];
                ++a; ++b; 
             }
             sscanf(dummystring, "%s", ligand_coords[i].atomname);

	for(a=0;a<20;++a)
                dummystring[a] = '\0';
             a=0; b = 17;
             while(b<=19)
             {   
                dummystring[a] =  line[b];
                ++a; ++b; 
             }
             sscanf(dummystring, "%s", ligand_coords[i].residuetype);

		for(a=0;a<20;++a)
                dummystring[a] = '\0';
             a=0; b = 22;
             while(b<=25)
             {   
                dummystring[a] =  line[b];
                ++a; ++b; 
             }
             sscanf(dummystring, "%d", &ligand_coords[i].seq_position);

		
		for(a=0;a<20;++a)
                dummystring[a] = '\0';
             a=0; b = 30;
             while(b<=37)
             {   
                dummystring[a] =  line[b];
                ++a; ++b; 
             }
                sscanf(dummystring, "%lf", &ligand_coords[i].coord.x);
             
             for(a=0;a<20;++a)
                dummystring[a] = '\0';
             a=0; b = 38;
             while(b<=45)
             {   
                dummystring[a] =  line[b];
                ++a; ++b; 
             }
                sscanf(dummystring, "%lf", &ligand_coords[i].coord.y);       

             for(a=0;a<20;++a)
                dummystring[a] = '\0';
             a=0; b = 46;
             while(b<=53)
             {   
                dummystring[a] =  line[b];
                ++a; ++b; 
             }
                sscanf(dummystring, "%lf", &ligand_coords[i].coord.z);

		ligand_coords[i].sasa=0; 
		ligand_coords[i].born_radius=0;
		ligand_coords[i].atom_ptr=dummy_atom_ptr;
		
		++i;
		

	   }
	}
	

  strcpy(ligand_coords[i].residuetype,"END");
  ligand_coords[i].sasa=ENDFLAG;
  ligand_coords[i].seq_position=ENDFLAG;
  ligand_coords[i].atom_number=ENDFLAG;
  ligand_coords[i].atom_ptr = NULL;


	free_memory(line);
	free_memory(dummystring);
}


/* read the ligand data in file ligandfilename. Place ligand info into the main resparam and rotamerlib structures.
	This function is called from read_forcefield */

void read_ligand_data(char *ligandfilename, RESPARAM *resparam, ROTAMERLIB *rotamerlib,  ATOMPARAM *atomparam)
{
	extern int GET_PID;
	FILE *ligand_file_ptr, *dummy_file_ptr;
	char *dummyfilename, *line, *dummystring;
	RESPARAM *working_resparam;
	ROTAMERLIB *working_rotamerlib;
	int i,j,k,u;
	extern double *STATIONARY_LIGAND_TRANSFORM_MATRIX;
	

	/* STATIONARY_LIGAND_TRANSFORM_MATRIX used when a dummy set of "dihedrals" are needed; ie: born_bonded */
	STATIONARY_LIGAND_TRANSFORM_MATRIX = (double *)calloc(14,sizeof(double));
	STATIONARY_LIGAND_TRANSFORM_MATRIX[1] = 1;	/* x */
	STATIONARY_LIGAND_TRANSFORM_MATRIX[2] = 0;
	STATIONARY_LIGAND_TRANSFORM_MATRIX[3] = 0;
	STATIONARY_LIGAND_TRANSFORM_MATRIX[4] = 0;

	STATIONARY_LIGAND_TRANSFORM_MATRIX[5] = 0;
	STATIONARY_LIGAND_TRANSFORM_MATRIX[6] = 1;	/* y */
	STATIONARY_LIGAND_TRANSFORM_MATRIX[7] = 0;
	STATIONARY_LIGAND_TRANSFORM_MATRIX[8] = 0;

	STATIONARY_LIGAND_TRANSFORM_MATRIX[9] = 0;
	STATIONARY_LIGAND_TRANSFORM_MATRIX[10] = 0;
	STATIONARY_LIGAND_TRANSFORM_MATRIX[11] = 1;	/* z */
	STATIONARY_LIGAND_TRANSFORM_MATRIX[12] = 0;

	working_resparam = (RESPARAM *)calloc(MAXRESTYPES,sizeof(RESPARAM));
	working_rotamerlib = (ROTAMERLIB *)calloc(MAXRESTYPES,sizeof(ROTAMERLIB));

	line = (char *)calloc(MXLINE_INPUT,sizeof(char));
	dummystring = (char *)calloc(MXLINE_INPUT,sizeof(char));

	dummyfilename = (char *)calloc(MXLINE_INPUT,sizeof(char));
	sprintf(dummyfilename,"%d.resrot",GET_PID);

	ligand_file_ptr = fopen_file(ligandfilename,"r");

	/* read resparam for ligand */

	fgets(line,MXLINE_INPUT,ligand_file_ptr);			/* move past notes */
	while(strncmp(line,"START",5)!=0)
		fgets(line,MXLINE_INPUT,ligand_file_ptr);
	
	dummy_file_ptr = fopen_file(dummyfilename,"w");		/* copy resparam entry to dummyfile */
	while(strncmp(line,"STOP",4)!=0)
	{
		fprintf(dummy_file_ptr,"%s",line);
		fgets(line,MXLINE_INPUT,ligand_file_ptr);
	}
	fprintf(dummy_file_ptr,"%s",line);
	fclose(dummy_file_ptr);

	read_residuedata(dummyfilename, working_resparam, atomparam, 0);	/* read resparam from dummyfile */

	working_resparam[1].ligand_flag=1;

	u=1;	
	while(strcmp(atomparam[u].atomtype,"U")!=0)
		++u;

	/* see whether bkbn atoms are defined; if not, give error */
	/* need to eventually add UUU bkbn as neeeded, and renumber the atomnumbers appropriately */

	if(strcmp(working_resparam[1].atom[4].atomname,"N")!=0)	/* backbone not defined! */
	{
		fprintf(stderr,"ERROR Need to define the resparam entry for %s with virtual backbone atoms!\n",working_resparam[1].residuetype);
		exit(1);
	}
	else	/* set backbone atoms to type U */
	{
		k=4;
		while(strcmp(working_resparam[1].atom[k].atomname,"zzz")!=0)
		{
			if(strcmp(working_resparam[1].atom[k].atomname,"N")==0 ||
				strcmp(working_resparam[1].atom[k].atomname,"CA")==0 ||
				strcmp(working_resparam[1].atom[k].atomname,"CB")==0 ||
				strcmp(working_resparam[1].atom[k].atomname,"O")==0 || strcmp(working_resparam[1].atom[k].atomname,"C")==0 ||
				strcmp(working_resparam[1].atom[k].atomname,"H")==0 || strcmp(working_resparam[1].atom[k].atomname,"HA")==0 ||
				working_resparam[1].atom[k].atomname[0] == 'U')
				{
					strcpy(working_resparam[1].atom[k].atomtype,"U");
					working_resparam[1].atom[k].atom_ptr = &atomparam[u];
				} 
			++k;
		}
	}

	
	/* read "rotamer" library */ 

	fgets(line,MXLINE_INPUT,ligand_file_ptr);			/* move past notes */
	while(strncmp(line,"START",5)!=0)
		fgets(line,MXLINE_INPUT,ligand_file_ptr);
	
	dummy_file_ptr = fopen_file(dummyfilename,"w");		/* "rotamer" entry to dummyfile */
	while(strncmp(line,"STOP",4)!=0)					
	{
		fprintf(dummy_file_ptr,"%s",line);
		fgets(line,MXLINE_INPUT,ligand_file_ptr);
	}
	fprintf(dummy_file_ptr,"%s",line);
	fclose(dummy_file_ptr);
	read_rotamers(dummyfilename,working_rotamerlib);

	/* read ligand base coordinates; superceeded by ligand coords in the Template file (if any) */
	if(fgets(line,MXLINE_INPUT,ligand_file_ptr)!=NULL)			/* move past notes */
	{
		while(strncmp(line,"START",5)!=0)
			fgets(line,MXLINE_INPUT,ligand_file_ptr);
		
		working_rotamerlib[1].ligand_base_coords = (pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM));
		read_ligand_coords(ligand_file_ptr,working_rotamerlib[1].ligand_base_coords);  /* read ligand_base_coords */
	}
	
	fclose(ligand_file_ptr);

	if(working_rotamerlib[1].ligand_base_coords==NULL)
	{
		if(working_rotamerlib[1].rotamer[1].ligand_coords==NULL)
		{
			fprintf(stderr,"ERROR Coordinates must be defined for ligand %s either as 'rotamers' or as base coordinates\n",
						working_rotamerlib[1].residuetype);
			exit(1);
		}
	}
	else
	{	
		if(working_rotamerlib[1].rotamer[1].ligand_coords!=NULL)	/* copy base coords to rotamerlib */
		{
			++working_rotamerlib[1].number_of_rotamers;
			working_rotamerlib[1].rotamer[working_rotamerlib[1].number_of_rotamers].ligand_coords = 
					(pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM));

			i=1;
			while(working_rotamerlib[1].ligand_base_coords[i].seq_position!=ENDFLAG)
			{
				working_rotamerlib[1].rotamer[working_rotamerlib[1].number_of_rotamers].ligand_coords[i] = 
								working_rotamerlib[1].ligand_base_coords[i];
				++i;
			}
			
			free_memory(working_rotamerlib[1].ligand_base_coords);
		}
	}
	working_resparam[1].rotamerlib_ptr = &working_rotamerlib[1];
	

	/* attach ligand resparam and rotamers to the main resparam and rotamer arrays */

	j=1;
	while(strcmp(resparam[j].residuetype,"zzz")!=0)
		++j;

	k=1;
	while(strcmp(rotamerlib[k].residuetype,"zzz")!=0)
		++k;

	resparam[j] = working_resparam[1];
	rotamerlib[k] = working_rotamerlib[1];
		
	rm_file(dummyfilename);

	free_memory(working_resparam);
	free_memory(working_rotamerlib);
	
	free_memory(dummyfilename);
	free_memory(line);
	free_memory(dummystring);
}



/* This function is called after the Template pdb file has been read, and the ligandparam files loaded,
	but before the VARIABLE_POSITIONS have been parsed.
	It creates virtual backbone atoms for each ligand (for book-keeping purposes).
	These (as well as the ligand atoms) are considered to be in chain 'l'
	It extracts the ligand coords from the Template pdb file (if any)
		if "rotamers" = transform matrix, coordinate from pdb file = base_coords for transformation
	Calculates intrinsic energies for ligands and protein sidechains */	

void initialize_stuff_for_ligands(PROTEIN *protein)
{
	char *ligand_tempfilename;
	CARTESIAN centroid, rotated_centroid, delta_centroid_vector;
	double *rotation_matrix;
	FILE *ligand_file_ptr;
	int i,r,n,k, template_pdb_ctr, seq_position, superchain_ctr;
	extern int GET_PID;
	pdbATOM *bkbn_pdb, *ligand_pdb;
	ENERGY free_ligand_energy;
	pdbATOM *ligand_res_pdb;
	mini_pdbATOM *ligand_res_minipdb;
	extern int MAX_SEQ_POSITION;
	div_t chain_id;

	ligand_tempfilename = (char *)calloc(MXLINE_INPUT, sizeof(char));
	bkbn_pdb = (pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM));
	ligand_res_minipdb = (mini_pdbATOM *)calloc(MAX_RES_SIZE,sizeof(mini_pdbATOM));
	ligand_res_pdb = (pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM));

	if(protein->chain_id_list == NULL)
	{
			protein->chain_id_list = (char *)calloc(MAX_CHAINS, sizeof(char));   
			i=0;
			while(i<MAX_CHAINS)
			{
				protein->chain_id_list[i] = ' '; 
				++i;
			}
			protein->chain_id_list[MAX_CHAINS] = '\0';
	 }

	if(protein->super_chain_id_list == NULL)
	{
		protein->super_chain_id_list = (SUPER_CHAIN_ID_LIST *)calloc(MAX_CHAINS,sizeof(SUPER_CHAIN_ID_LIST));
                        for(i=1;i<MAX_CHAINS;++i)
                                protein->super_chain_id_list[i].chain_id=NULL;

                protein->super_chain_id_list[1].chain_id = (char *)calloc(MAX_CHAINS,sizeof(char));
                protein->super_chain_id_list[1].chain_id[1] = protein->chain_id_list[1];
                protein->super_chain_id_list[1].chain_id[2] = '\0';
	}

	seq_position = protein->seqpos_text_map[protein->num_res].seq_position;

			chain_id = div(seq_position, 1000);
			seq_position = 1000*(chain_id.quot + 1)+1; /* start ligands as new chain l */
			chain_id = div(seq_position, 1000);
			protein->chain_id_list[chain_id.quot+1] = 'l';
			protein->chain_id_list[chain_id.quot+2] = '\0';

	superchain_ctr=1;
	while(protein->super_chain_id_list[superchain_ctr].chain_id!=NULL)
		++superchain_ctr;		
	protein->super_chain_id_list[superchain_ctr].chain_id = (char *)calloc(MAX_CHAINS,sizeof(char));
        protein->super_chain_id_list[superchain_ctr].chain_id[1] = protein->chain_id_list[chain_id.quot+1];
        protein->super_chain_id_list[superchain_ctr].chain_id[2] = '\0';

	
	r=1; template_pdb_ctr=1; 
	while(strcmp(protein->resparam[r].residuetype,"zzz")!=0)
	{
		if(protein->resparam[r].ligand_flag==1)
		{
			sprintf(ligand_tempfilename,"temp.%d.%s.lig",GET_PID,protein->resparam[r].residuetype);

			if(does_this_file_exist(ligand_tempfilename)==1)	/* read in coords from Template pdb file */
			{
				ligand_file_ptr = fopen_file(ligand_tempfilename,"r");
				
			/* since ligand_base_coords!=NULL, using transform matrix */
				if(protein->resparam[r].rotamerlib_ptr->ligand_base_coords!=NULL)	/* replace base coords */
				{
					read_ligand_coords(ligand_file_ptr,protein->resparam[r].rotamerlib_ptr->ligand_base_coords); 
	
					if(protein->transform_matrix!=NULL)
						transform_pdb(protein->resparam[r].rotamerlib_ptr->ligand_base_coords, protein->transform_matrix);

				}
				else  /* since the base coords are NULL, the actual coords are being used as rotamers */
				{
					++protein->resparam[r].rotamerlib_ptr->number_of_rotamers;

		protein->resparam[r].rotamerlib_ptr->rotamer[protein->resparam[r].rotamerlib_ptr->number_of_rotamers].ligand_coords = 
						(pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM));

				read_ligand_coords(ligand_file_ptr,protein->resparam[r].rotamerlib_ptr->rotamer[protein->resparam[r].rotamerlib_ptr->number_of_rotamers].ligand_coords); 


				if(protein->transform_matrix!=NULL)
					transform_pdb(protein->resparam[r].rotamerlib_ptr->rotamer[protein->resparam[r].rotamerlib_ptr->number_of_rotamers].ligand_coords, protein->transform_matrix);

					
				}					

				fclose(ligand_file_ptr); 

				rm_file(ligand_tempfilename);
			}

			// ligamers were automatically generated, but the rotation vectors are based at the origin; 
			// transfer them to the centroid of the base coordinates
			if(protein->resparam[r].rotamerlib_ptr->generate_ligamers_flag==1)
			{
				rotation_matrix = (double *)calloc(14,sizeof(double));	

				CARTESIAN_zero(centroid);
        			i=0; k=1;
        			while(protein->resparam[r].rotamerlib_ptr->ligand_base_coords[k].seq_position!=ENDFLAG)
        			{
                			if(strcmp(protein->resparam[r].rotamerlib_ptr->ligand_base_coords[k].atomname,"U")!=0)
                			{
                        		centroid = 
							add_CARTESIAN(centroid, protein->resparam[r].rotamerlib_ptr->ligand_base_coords[k].coord);
	
    	                  		++i;
                			}
                			++k;
        			}

        			centroid = scalar_multiply_CARTESIAN((double)(1.0/i), centroid);
				
				/*
				printf("%lf\t%lf\t%lf\n",centroid.x,centroid.y,centroid.z);
			 
				for(i=1;i<=12;++i)
					printf("%.18lf ",protein->resparam[r].rotamerlib_ptr->rotamer[2].chi[i]);
				printf("\n\n");
			    */

				// adjust transform matricies 
			
				for(i=1;i<=protein->resparam[r].rotamerlib_ptr->number_of_rotamers;++i)
				{
					for(k=1;k<=12;++k)
						rotation_matrix[k] = protein->resparam[r].rotamerlib_ptr->rotamer[i].chi[k];
					rotation_matrix[4] = 0;
					rotation_matrix[8] = 0;
					rotation_matrix[12] = 0;

					rotated_centroid = tranform_coordinates(centroid, rotation_matrix);
					
					delta_centroid_vector = subtract_CARTESIAN(centroid, rotated_centroid);
				

					protein->resparam[r].rotamerlib_ptr->rotamer[i].chi[4] += delta_centroid_vector.x;
					protein->resparam[r].rotamerlib_ptr->rotamer[i].chi[8] += delta_centroid_vector.y;
					protein->resparam[r].rotamerlib_ptr->rotamer[i].chi[12] += delta_centroid_vector.z;
				}

			 	 /*
				for(i=1;i<=12;++i)
					printf("%.18lf ",protein->resparam[r].rotamerlib_ptr->rotamer[2].chi[i]);
				printf("\n");
				*/
			    

				sprintf(ligand_tempfilename,"ligamers.%s.%d",protein->resparam[r].rotamerlib_ptr->residuetype,GET_PID);
				rm_file(ligand_tempfilename);

				free_memory(rotation_matrix);

			}


			/* generate backbone atoms */			

			chain_id = div(seq_position, 1000);
			++protein->num_res;

			sprintf(protein->seqpos_text_map[protein->num_res].seqpos_text,"%dl",chain_id.rem);
			protein->seqpos_text_map[protein->num_res].seq_position = seq_position;
		
			if(protein->seqpos_text_map[protein->num_res].seq_position > MAX_SEQ_POSITION)
  				MAX_SEQ_POSITION = protein->seqpos_text_map[protein->num_res].seq_position;
  			protein->seqpos_text_map[protein->num_res+1].seq_position = ENDFLAG;
			
			
			while(strcmp(protein->Template[template_pdb_ctr].residuetype,"PRO")==0 || 
					strcmp(protein->Template[template_pdb_ctr].residuetype,"GLY")==0 ||
					strcmp(protein->Template[template_pdb_ctr].residuetype,"CYD")==0)
				++template_pdb_ctr;

			n=1; k=protein->Template[template_pdb_ctr].seq_position;
			while(protein->Template[template_pdb_ctr].seq_position == k)
			{
				if(protein->Template[template_pdb_ctr].atom_ptr->other_info<0)
				{
					bkbn_pdb[n] = protein->Template[template_pdb_ctr];
					++n;
				}
				++template_pdb_ctr;
			}
			bkbn_pdb[n].atom_ptr=NULL; 	bkbn_pdb[n].seq_position=ENDFLAG; 	
			strcpy(bkbn_pdb[n].residuetype,"END");	
			
			for(n=1;n<=7;++n)
			{
	    			bkbn_pdb[n].atom_ptr = &protein->resparam[r].sorted_atom[search_atomname_ATOMRESPARAM(bkbn_pdb[n].atomname,
																				protein->resparam[r].sorted_atom)];
	    			bkbn_pdb[n].seq_position = seq_position;
	    			bkbn_pdb[n].born_radius=0;
	    			bkbn_pdb[n].sasa = 0;
				bkbn_pdb[n].atom_number = n;
				strcpy(bkbn_pdb[n].residuetype,protein->resparam[r].residuetype);
				strcpy(bkbn_pdb[n].seqpos_text,protein->seqpos_text_map[protein->num_res].seqpos_text);
	    		}
	    	
			/* put correct seq_position for ligand_coord or ligand_base_coord */
			if(protein->resparam[r].rotamerlib_ptr->ligand_base_coords!=NULL)
			{
				n=1;
				while(protein->resparam[r].rotamerlib_ptr->ligand_base_coords[n].seq_position!=ENDFLAG)
				{
					protein->resparam[r].rotamerlib_ptr->ligand_base_coords[n].seq_position = seq_position;
					strcpy(protein->resparam[r].rotamerlib_ptr->ligand_base_coords[n].seqpos_text,
										protein->seqpos_text_map[protein->num_res].seqpos_text);
					++n;			
				}
				ligand_pdb = protein->resparam[r].rotamerlib_ptr->ligand_base_coords;
			}
			else
			{
				for(i=1;i<=protein->resparam[r].rotamerlib_ptr->number_of_rotamers;++i)
				{
					n=1;
					while(protein->resparam[r].rotamerlib_ptr->rotamer[i].ligand_coords[n].seq_position!=ENDFLAG)
					{
						protein->resparam[r].rotamerlib_ptr->rotamer[i].ligand_coords[n].seq_position = seq_position;
						++n;			
					}
				}

				ligand_pdb = 
			protein->resparam[r].rotamerlib_ptr->rotamer[protein->resparam[r].rotamerlib_ptr->number_of_rotamers].ligand_coords;
				--protein->resparam[r].rotamerlib_ptr->number_of_rotamers; 
				/* since the last rotamer is being attached to the template, it will be re-extracted later */
			}

			/* append backbone to Template pdb */
			append_pdbATOM(bkbn_pdb,protein->Template);

			/* append ligand atoms to the Template pdb */  
			append_pdbATOM(ligand_pdb,protein->Template);

			/* create ligand res and evaluate free ligand energy */
			ligand_res_pdb[1].atom_ptr=NULL; ligand_res_pdb[1].seq_position=ENDFLAG; 
			strcpy(ligand_res_pdb[1].residuetype,"END");
			append_pdbATOM(bkbn_pdb,ligand_res_pdb);
			append_pdbATOM(ligand_pdb,ligand_res_pdb);

			attach_ptr_pdbATOM(ligand_res_pdb, protein->resparam);


			pdbATOM_to_mini_pdbATOM(ligand_res_pdb, ligand_res_minipdb);
			free_ligand_energy = energy_calc(ligand_res_minipdb,1,1);
			if(protein->resparam[r].E_random == 0)
			{
				protein->resparam[r].E_random = ENERGY_TOTAL(free_ligand_energy);
				protein->resparam[r].E_unfolded = protein->resparam[r].E_random;
				protein->resparam[r].E_reference = protein->resparam[r].E_random;
				protein->resparam[r].E_tripeptide = protein->resparam[r].E_random;
			}	
			
			++seq_position; 

		}
		++r;
	}

	 
	attach_ptr_pdbATOM(protein->Template, protein->resparam);

	 
	pdbATOM_to_mini_pdbATOM(protein->Template, protein->mini_Template); 

	i=1; protein->num_res=0;
        while(protein->mini_Template[i].seq_position!=ENDFLAG)
        {
                if(strcmp(protein->mini_Template[i].atom_ptr->atomname,"N")==0)
                        ++protein->num_res;
                ++i;
        }

	free_memory(ligand_tempfilename); free_memory(bkbn_pdb); free_memory(ligand_res_minipdb); free_memory(ligand_res_pdb);

	intrinsic_rotamer_energy(protein->resparam);	
}