/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
**                                                                                 **
**                                                                                 **
**                                                       **                        **
**                                                       **                        **
**                       *********                       **                        **
**        **********    **********    ***        ***     **                        **
**        **********    **            ****      ****     **      **********        **
**        **            **            **  **  **  **     **      **********        **
**        **            ********      **  **  **  **     **      **                **
**        ******         ********     **    **    **     **      **                **
**        ******               ***    **    **    **     **************            **
**        **     ***     **    **********         **     **************            **
**        **     ***    *****************        ********        **                **
**        ***********   *********     **         *********       **                **
**        ***********    **    **     **         ***    **       **********        **
**               **  **  **    ******            **     **       **********        **
**               **  **  **    ******            ********                          **
**               **    * **    **                *********                         **
**               **    * **    **                **      **                        **
**               **     ***    **********        **      **                        **
**               **     ***    **********        **********                        **
**                                               *********                         **
**                                                                                 **
**                                                                                 **
**                          E N S E M B L E     T O O L S                          **
**                                                                                 **
**                                   Version 2.1                                   **
**                                                                                 **
**                                                                                 **
**                   Principal Investigator: Julie D. Forman-Kay                   **
**                                                                                 **
**                           Author:  Mickaël Krzeminski                           **
**                                                                                 **
**                               Date: November 2012                               **
**                                                                                 **
**                                                                                 **
*************************************************************************************
**                                                                                 **
**                                                                                 **
**   Copyright (C) The Hospital for Sick Children, 2001                            **
**                                                                                 **
**   Distribution of substantively modified versions of this module is prohibited  **
**  without the explicit permission of the copyright holder.                       **
**                                                                                 **
**   Any use of this work or derivative works in whole or in part for any          **
**  commercial purpose or for monetary gain is prohibited.                         **
**                                                                                 **
**                                                                                 **
**                                   NO WARRANTY                                   **
**   This software package is provided 'as is' without warranty of any kind,       **
**  expressed or implied.                                                          **
**                                                                                 **
**                                                                                 **
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef __ACCESSURF_HEADER__
#define __ACCESSURF_HEADER__


counter errors;

void _cleanAS(PDBListData *ENSData, accessSurf *access, char *PDBFileName, char *outFileName, char *FFName)
{
	// This functions releases the memory affected in the main function
	_cleanPDB(ENSData);

	if (access != NULL)
	{
		free((*access).as);
		(*access).as = NULL;
	}

	free(PDBFileName);
	free(outFileName);
	free(FFName);
}


void _cleanh(int ict, unsigned int *tag, unsigned int *tag1, int *kn, int *inz, double *zlb, double *zub, double *rad, double *rad1, double *xr1, double *yr1, double *zr1,\
            double *rsec, double *rsec2, double *area, double **arci, double **arcf, protein pdb, typerad *radia)
{
	// This functions releases the memory affected in the heading function getAccSurf
	counter count;

	free(tag); free(tag1);
	free(kn);
	free(inz);
	free(zlb); free(zub);
	free(rad); free(rad1);
	free(xr1); free(yr1); free(zr1);
	free(area); //free(access);
	free(rsec);
	free(rsec2);
	free(radia);

	for (count=0; count<pdb.nbResidues; count++)
	{
		free(pdb.resAtoms[count].names);
		free(pdb.resAtoms[count].types);
	}
	free(pdb.resAtoms);

	if (arci != NULL)
	{
		for (count=0; count<ict; count++) free(arci[count]);
	}
	if (arcf != NULL)
	{
		for (count=0; count<ict; count++) free(arcf[count]);
	}

	free(arci); free(arcf);

	return;
}


void _showError(char *var)
{
	char msg[40];
	printf ("VAR = %s\n", var);
	strcpy(msg, "instance ");
	strcpy(&msg[9], var);
	strcpy(&msg[9+strlen(var)], " variable");
	error("MEMORY", msg);
	++errors;
}


// The following sortig table is based on the quicksort algorithm from http://interstices.info/jcms/c_6973/les-algorithmes-de-tri
void swapdouble(double *values, int a, int b)
{
	double temp = values[a];
	
	values[a] = values[b];
	values[b] = temp;
}


void swapuint32(uint32_t *values, int a, int b)
{
	uint32_t temp = values[a];
	
	values[a] = values[b];
	values[b] = temp;
}


unsigned int partition(double *values, uint32_t *classif, const int beg, const int end)
{
	counter i, count = beg;
	uint32_t pivotIndex;

	srand((unsigned int)time(NULL));
	pivotIndex = rand()%(end-beg+1)+beg;		// Random values within [beg, end]
	swapdouble(values, beg, pivotIndex);
	swapuint32(classif, beg, pivotIndex);

	for (i=beg+1; i<=end; i++)
	{
		if (values[i] < values[beg])
		{
			++count;
			swapdouble(values, count, i);
			swapuint32(classif, count, i);
		}
	}

	swapdouble(values, count, beg);
	swapuint32(classif, count, beg);
	return count;
}


void fastSortAux(double *values, uint32_t *classif, const int beg, const int end)
{
	if (beg < end)
	{
		unsigned int pivotIndex = partition(values, classif, beg, end);
		fastSortAux(values, classif, beg, pivotIndex-1);
		fastSortAux(values, classif, pivotIndex+1, end);
	}
}


void sortag(double *values, uint32_t *classif, unsigned int len)
{
	counter i;

	for (i=0; i<len; i++) classif[i] = i;

	fastSortAux(values, classif, 0, len-1);
}


void calarc(int m, int isign, double d2, double q, double d, double *rsec, double **arci, double **arcf, int *kn, double dy, double dx)
{
	int k1;
	double arg, alpha1, beta1, ti, tf;
	double piex2 = 2 * M_PI;


	k1 = kn[m];
	// Law of cosines
	// alpha1 is the angle between a line containing a point of intersection and
	// the reference circle center and the line containing both circle centers
	arg = (d2 + q * isign) / (d * rsec[m]);
	alpha1 = acos(arg);

	// beta1 is the angle between the line containing both circle centers and the x-axis
	beta1 = atan2(dy, dx);

	if (isign == 1) beta1 += M_PI;
	ti = beta1 - alpha1;
	tf = beta1 + alpha1;

	if (ti < 0.) ti += piex2;
	if (tf > piex2) tf -= piex2;
	if (tf < 0.) tf += piex2;

	if (tf < ti)
	{
		arcf[k1+1][m] = tf;
		arcf[k1][m] = piex2;
		++kn[m];
	}
	else
	{
		arcf[k1][m] = tf;
	}

	arci[k1][m] = ti;
}


void getColumn(double **origin, double *desti, unsigned int column, unsigned int length)
{
	counter i;

	for (i=0; i<length; i++) desti[i] = origin[i][column];
}


void setColumn(double *origin, double **desti, unsigned int column, unsigned int length)
{
	counter i;

	for (i=0; i<length; i++) desti[i][column] = origin[i];
}


int reduce(int *n, int itab, int ict, double **arci, double **arcf, uint32_t *tag, uint32_t *tag1, double *zr1, double *rad1, double *area, int *inz,\
           int *kn, int *icnt1, int *icnt2, int ict1, double zres, double rh2o, double zgrid, double hzres, double cutoff)
{
	int ii, jj, k, m, i, i1, j, k11, k1 = 0;
	double arcsum, t, a, parea = 0.;
	double piex2 = 2. * M_PI;
	double *colArci = malloc(DOUBLE_SIZE);

	for (ii=0; ii<max(1, itab); ii++)
	{
		if (itab != 0) (*n) = inz[ii];
		k1 = kn[(*n)];
		if (k1 > ict-1) continue;

		if (k1 == -1)
		{
			arcsum = piex2;
		}
		else
		{
			// The arc endpoints are sorted on the value of the initial arc endpoint
			k11 = k1 + 1;
			if ((colArci = realloc(colArci, k11 * DOUBLE_SIZE)) == NULL) return 1;

			getColumn(arci, colArci, (*n), k11);
			sortag(colArci, tag, k11);
			setColumn(colArci, arci, (*n), k11);
			
			// Calculate the accessible area
			arcsum =  arci[0][(*n)];
			jj = tag[0];
			t = arcf[jj][(*n)];
			if (k1 != 0)
			{
				for (k=1; k<max(2, k11); k++)
				{
					if (t < arci[k][(*n)])
					{
						arcsum += arci[k][(*n)] - t;
						jj = tag[k];
						t = arcf[jj][(*n)];
					}
					else
					{
						m = tag[k];
						if (arcf[m][(*n)] > t)
						{
							t = arcf[m][(*n)];
						}
					}
				}
			}
			arcsum += piex2 - t;
		}
		a = zr1[(*n)] - zgrid;
		a = rad1[(*n)] - fabs(a);

		// The area is equal to the accessible arc length times the section thickness,
		// corrected if it is the first or last section, times the radius of the sphere
		parea = arcsum * (hzres + min(a, hzres)) * rad1[(*n)];

		// Add the accessible area for this atom in this section to the area
		// for this atom for all the section encountered thus far
		area[(*n)] += parea;

	}
	free(colArci);

	// If this was the final atom for the z section return
	if (itab > 0) return 0;

	// The arci and arcf arrays were filled, does significant area remain ?
	if (parea <= cutoff)
	{
		kn[(*n)] = 9999;
		++(*icnt1);
		return 0;
	}

	// Removing degeneracies from arci and arcf (Combine overlapping arcs)
	jj = tag[0];
	t = arcf[jj][(*n)];
	i = 0;

	for (k=1; k<max(2, k1+1); k++)
	{
		if (t < arci[k][(*n)])
		{
			++i;
			if (i == ict1-1)
			{
				printf("Error in surface calculation: Atom %d has more than %d contacts!\n\n", tag1[(*n)], ict);
				return 1;
			}

			i1 = i - 1;
			for (j=k; j<max(k+1, k1+1); j++)
			{
				if (tag[j] == i1)
				{
					arcf[jj][(*n)] = arcf[i1][(*n)];
					tag[j] = jj;
					break;
				}
			}
			arcf[i1][(*n)] = t;
			arci[i][(*n)] = arci[k][(*n)];
			jj = tag[k];
			t = arcf[jj][(*n)];
		}
		else
		{
			m = tag[k];
			if (arcf[m][(*n)] > t)
			{
				t = arcf[m][(*n)];
			}
		}
	}

	arcf[i][(*n)] = t;
	kn[(*n)] = i;
	++i;
	for (k=i; k<max(i+1, k1+1); k++) arci[k][(*n)] = 0.;

	// Removing the partial area added, this circle may have more intersections
	area[(*n)] -= parea;
	++(*icnt2);
	return 0;
}


accessSurf getAccSurf(PDBListData *ENSData, char *FFName, double rh2o, char mode, char radtype, double accuracy)
{
	// Variable declaration
	extern counter errors;
	counter i, j, k;
	counter nbTypeRad = 0;
	counter maxNbResidues = 50, maxNbAtoms = 0;

	enum {False, True} resok = False, endok = True, ffok = False;
	const int ict = 15;
	const unsigned int dim = max(ict, (*ENSData).nbAtoms);

	char uc;
	char atypeFileLine[ATYPE_LINE_LEN];
	char valchar[20], atypeElement[8], ff[20];

	unsigned int ucseek, ffCol = 1;
	unsigned int *tag;
	unsigned int *tag1;

	int ict1 = ict - 1;
	int termres = -1;
	int iano, j1, iflag;
	int itab, l, m, idum=0;
	int *kn;
	int *inz;

	double zubmax = 0., hzres, zor, znext;
	double dx, dy, d, d2, q, b;
	double rmin = 10000., cutoff, zgrid, a;
	double radFact= 1./pow(2, 5./6.);

	double *zlb, *zub;		// Arrays that contain z coordinates of upper and lower bounds
	double *rad,  *rad1;		// Arrays that contain the radius of all atoms
	double *xr1, *yr1, *zr1;	// Arrays that contain translated x, y, and z coordinates
	double *rsec, *rsec2;
	double *area;

	double **arci;
	double **arcf;


	protein pdb;
	typerad *radia;
	accessSurf access;

	FILE *atypeFile;

	// Array variable assignement
	errors = 0;
	access.FAILED = 1;

	if ((tag       = malloc(dim*INT_SIZE))    == NULL) _showError("tag");
	if ((tag1      = malloc(dim*INT_SIZE))    == NULL) _showError("tag1");
	if ((kn        = malloc(dim*INT_SIZE))    == NULL) _showError("kn");
	if ((inz       = malloc(dim*INT_SIZE))    == NULL) _showError("inz");
	if ((zlb       = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("zlb");
	if ((zub       = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("zub");
	if ((rad       = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("rad");
	if ((rad1      = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("rad1");
	if ((xr1       = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("xr1");
	if ((yr1       = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("yr1");
	if ((zr1       = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("zr1");
	if ((rsec      = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("rsec");
	if ((rsec2     = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("rsec2");
	if ((area      = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("area");
	if ((access.as = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("access.as");

	if ((arci     = malloc(ict*POINTER_SIZE)) == NULL) _showError("arci");
	if ((arcf     = malloc(ict*POINTER_SIZE)) == NULL) _showError("arcf");
	if (arci != NULL && arcf != NULL)
	{
		for (i=0; i<ict; i++)
		{
			if ((arci[i] = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("arci[]");
			if ((arcf[i] = malloc(dim*DOUBLE_SIZE)) == NULL) _showError("arcf[]");
		}
	}

	if ((pdb.resAtoms = malloc(maxNbResidues*sizeof(residue))) == NULL) _showError("pdb.resAtoms");

	if ((radia = malloc(100*sizeof(typerad))) == NULL) _showError("radia");

	if (errors)
	{
		_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
		return access;
	}


	// Recording the atoms radia, as well as the zlb array values
	if ((atypeFile = fopen("ensatom.types", "rt")) == NULL)
	{
		error("NO_FILE", "ensatom.types");
		_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
		return access;
	}

	while (fgets(atypeFileLine, ATYPE_LINE_LEN, atypeFile) != NULL)
	{
		uc = 0;
		for (i=0; i<min(ATYPE_LINE_LEN, strlen(atypeFileLine)); i++)
		{
			uc = toupper(atypeFileLine[i]);
			if (uc>32 && uc<127) break;			// Any character found
		}

		if (uc<48 || (uc>57 && uc<65) || uc>91) continue;	// Ignoring empty lines or lines starting with non-letter character

		for (j=0; j<7 && ((uc>47 && uc<58) || (uc>64 && uc<91)); j++)
		{
			atypeElement[j] = uc;
			++i;
			uc = toupper(atypeFileLine[i]);
		}
		atypeElement[j] = '\0';
		
		if (! strcmp(atypeElement, "RESIDUE"))			// New residue found
		{
			if (endok == False)
			{
				error_msg("The keyword 'END' was not found!\n", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}
	
			if (pdb.nbResidues == maxNbResidues)
			{
				maxNbResidues += 20;
				if ((pdb.resAtoms = realloc(pdb.resAtoms, maxNbResidues*sizeof(residue))) == NULL)
				{
					_showError("pdb.residues");
					_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
					return access;
				}
			}

			while (uc<65 || uc>90)			// Reading until finding the residue name
			{
				++i;
				uc = toupper(atypeFileLine[i]);
			}

			for (j=0; j<3; j++)
			{
				pdb.resAtoms[pdb.nbResidues].resn[j] = uc;
				++i;
				uc = toupper(atypeFileLine[i]);
			}
			pdb.resAtoms[pdb.nbResidues].resn[j] = '\0';

			// Checking whether the residue name has not already been defined
			for (j=0; j<pdb.nbResidues; j++)
			{
				if (! strcmp(pdb.resAtoms[pdb.nbResidues].resn, pdb.resAtoms[j].resn))
				{
					error("Double definition of the same residue (%s)!\n\n", pdb.resAtoms[pdb.nbResidues].resn);
					_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
					return access;
				}
			}

			// Preparing the array for atom names and types
			pdb.resAtoms[pdb.nbResidues].nbAtoms = 0;
			maxNbAtoms = 50;
			if ((pdb.resAtoms[pdb.nbResidues].names = malloc(NTLEN*maxNbAtoms*CHAR_SIZE)) == NULL)
			{
				error("MEMORY", "instance pdb.resAtoms[pdb.nbResidues].names variable");
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}
			if ((pdb.resAtoms[pdb.nbResidues].types = malloc(NTLEN*maxNbAtoms*CHAR_SIZE)) == NULL)
			{
				error("MEMORY", "instance pdb.resAtoms[pdb.nbResidues].types variable");
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			resok = True;
			endok = False;
		}
		else if (! strcmp(atypeElement, "TERMINI"))			// Atom types that can be applied to all patched residues
		{
			if (endok == False)
			{
				error_msg("The keyword 'END' was not found!\n", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}
			
			if (pdb.nbResidues == maxNbResidues)
			{
				maxNbResidues += 20;
				if ((pdb.resAtoms = realloc(pdb.resAtoms, maxNbResidues*sizeof(residue))) == NULL)
				{
					fatal_error("MEMORY", "reinstance residues variable");
				}
			}

			// Preparing the array for atom names and types
			pdb.resAtoms[pdb.nbResidues].nbAtoms = 0;
			maxNbAtoms = 20;
			if ((pdb.resAtoms[pdb.nbResidues].names = malloc(NTLEN*maxNbAtoms*CHAR_SIZE)) == NULL)
			{
				_showError("pdb.resAtoms[pdb.nbResidues].names");
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}
			if ((pdb.resAtoms[pdb.nbResidues].types = malloc(NTLEN*maxNbAtoms*CHAR_SIZE)) == NULL)
			{
				_showError("pdb.resAtoms[pdb.nbResidues].types");
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			strcpy(pdb.resAtoms[pdb.nbResidues].resn, TERM_NAME);		// TERM_NAME is supposedly not used for any other residue
			termres = pdb.nbResidues;
			resok = True;
			endok = False;
		}
		else if (! strcmp(atypeElement, "END"))		// End of a residue
		{
			if (pdb.resAtoms[pdb.nbResidues].nbAtoms == 0)
			{
				error_msg("No atom specified for residue %s!\n", pdb.resAtoms[pdb.nbResidues].resn);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			if (resok == False)
			{
				error_msg("No residue specification before the END!\n", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			// Reallocation of pdb.resAtoms[pdb.nbResidues] to fit the number of atoms
			if ((pdb.resAtoms[pdb.nbResidues].names = realloc(pdb.resAtoms[pdb.nbResidues].names, NTLEN*pdb.resAtoms[pdb.nbResidues].nbAtoms*CHAR_SIZE)) == NULL)
			{
				error("MEMORY", "instance pdb.resAtoms[pdb.nbResidues].names variable");
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}
			if ((pdb.resAtoms[pdb.nbResidues].types = realloc(pdb.resAtoms[pdb.nbResidues].types, NTLEN*pdb.resAtoms[pdb.nbResidues].nbAtoms*CHAR_SIZE)) == NULL)
			{
				error("MEMORY", "instance pdb.resAtoms[pdb.nbResidues].types variable");
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			++pdb.nbResidues;
			endok = True;
			resok = False;
		}
		else						// New atom found
		{
			if (resok == False)
			{
				error_msg("No residue name found before atom enumeration!\n", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			if (pdb.resAtoms[pdb.nbResidues].nbAtoms == maxNbAtoms)
			{
				maxNbAtoms += 50;
				if ((pdb.resAtoms[pdb.nbResidues].names = realloc(pdb.resAtoms[pdb.nbResidues].names, NTLEN*maxNbAtoms*CHAR_SIZE)) == NULL)
				{
					error("MEMORY", "instance pdb.resAtoms[pdb.nbResidues].names variable");
					_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
					return access;
					
				}
				if ((pdb.resAtoms[pdb.nbResidues].types = realloc(pdb.resAtoms[pdb.nbResidues].types, NTLEN*maxNbAtoms*CHAR_SIZE)) == NULL)
				{
					error("MEMORY", "instance pdb.resAtoms[pdb.nbResidues].types variable");
					_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
					return access;
				}
			}

			if (strlen(atypeElement) > 4)
			{
				error("The name of an atom must not be longer that 4 characters!", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			strncpy(pdb.resAtoms[pdb.nbResidues].names[pdb.resAtoms[pdb.nbResidues].nbAtoms], atypeElement, 4);
			pdb.resAtoms[pdb.nbResidues].names[pdb.resAtoms[pdb.nbResidues].nbAtoms][4] = '\0';

			while (uc<48 || (uc>57 && uc<65) || uc>90)			// Reading until finding the atom type
			{
				++i;
				uc = toupper(atypeFileLine[i]);
			}

			for (j=0; j<4 && ((uc>47 && uc<58) || (uc>64 && uc<91) || uc==95); j++)
			{
				pdb.resAtoms[pdb.nbResidues].types[pdb.resAtoms[pdb.nbResidues].nbAtoms][j] = uc;
				++i;
				uc = toupper(atypeFileLine[i]);
			}
			pdb.resAtoms[pdb.nbResidues].types[pdb.resAtoms[pdb.nbResidues].nbAtoms][j] = '\0';
			
			if (uc>32 && uc<127)
			{
				error("The name of an atom type must not be longer that 4 characters!", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			++pdb.resAtoms[pdb.nbResidues].nbAtoms;
		}		
	}
	fclose(atypeFile);


	if ((atypeFile = fopen("ensigmatom.dat", "rt")) == NULL)
	{
		error("NO_FILE", "ensatom.types");
		_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
		return access;
	}

	while (fgets(atypeFileLine, ATYPE_LINE_LEN, atypeFile) != NULL)
	{
		uc = 0;
		for (i=0; i<min(ATYPE_LINE_LEN, strlen(atypeFileLine)); i++)
		{
			uc = toupper(atypeFileLine[i]);
			if (uc>32 && uc<127) break;			// Any character found
		}
		
		if (uc<48 || (uc>57 && uc<65) || uc>91) continue;	// Ignoring empty lines or lines starting with non-letter character

		for (j=0; j<4 && ((uc>47 && uc<58) || (uc>64 && uc<91) || uc==95); j++)
		{
			atypeElement[j] = uc;
			++i;
			uc = toupper(atypeFileLine[i]);
		}
		atypeElement[j] = '\0';

		if (! strcmp(atypeElement, "FFNA"))		// Forcefield list
		{
			if (ffok == True) continue;			// Ignoring other occurences of force field specification

			while (uc>64 && uc<91)				// Finishing reading the FFNAME keyword
			{
				++i;
				uc = toupper(atypeFileLine[i]);
			}

			while (1)					// Detection of the forcefield column
			{
				while (uc<33 || uc>126)				// Reading until finding the first character
				{
					++i;
					uc = toupper(atypeFileLine[i]);
					if (uc==10 || uc==13)
					{
						error("The specified forcefield was not found in the sigma value file!", NULL);
						_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
						return access;
					}
				}

				ucseek = 0;
				while (uc>32 && uc<127)				// Retrieving the forcefield name
				{
					ff[ucseek] = uc;
					++i; ++ucseek;
					uc = toupper(atypeFileLine[i]);
				}
				ff[ucseek] = '\0';

				if (! strcmp(ff, FFName))
				{
					ffok = True;
					break;
				}

				++ffCol;
			}
		}
		else
		{
			if (ffok == False)
			{
				error("No specified forcefield in the sigma value file!", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			strncpy(radia[nbTypeRad].type, atypeElement, 5);

			if (uc>32 && uc<127)
			{
				error("The name of an atom type must not be longer that 4 characters!", NULL);
				_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
				return access;
			}

			for (j=1; j<ffCol; j++)				// Going to the column of the selected forcefield
			{
				while (uc<46 || uc>57)			// Reading until finding a value
				{
					++i;
					uc = toupper(atypeFileLine[i]);
				}

				if (j==1)				// Recording the default value
				{
					ucseek = 0;
					while (uc>45 && uc<58)
					{
						valchar[ucseek] = atypeFileLine[i];
						++i; ++ucseek;
						uc = atypeFileLine[i];
					}
					valchar[ucseek] = '\0';

					if (valchar[0] != '/')
					{
						radia[nbTypeRad].radius = atof(valchar);
					}
					else
					{
						radia[nbTypeRad].radius = 0.;
					}
				}
				else
				{
					while (uc>45 && uc<58)
					{
						++i;
						uc = atypeFileLine[i];
					}
				}
			}

			while (uc<46 || uc>57)			// Reading until finding a value
			{
				++i;
				uc = toupper(atypeFileLine[i]);
			}


			ucseek = 0;
			while (uc>45 && uc<58)
			{
				valchar[ucseek] = atypeFileLine[i];
				++i; ++ucseek;
				uc = atypeFileLine[i];
			}
			valchar[ucseek] = '\0';

			if (valchar[0] != '/')
			{
				radia[nbTypeRad].radius = atof(valchar);
				++nbTypeRad;
			}
			else if (radia[nbTypeRad].radius > dEpsilon)
			{
				++nbTypeRad;
			}
		}
	}
	fclose(atypeFile);


	// Now attributing the radius to each atom of the molecule
	char name[5] = {'\0'};
	char resn[4] = {'\0'};
	char resi[5] = {'\0'};
	
	for (i=0; i<dim; i++)
	{
		rad[i] = 0.;
		strncpy(name, &(*ENSData).atomHeaders[i][12], 4);
		strncpy(resn, &(*ENSData).atomHeaders[i][17], 3);

		deepTrim(name); deepTrim(resn);

		for (j=0; j<pdb.nbResidues; j++)
		{
			if (! strcmp(pdb.resAtoms[j].resn, resn))
			{
				for (k=0; k<pdb.resAtoms[j].nbAtoms; k++)
				{
					if (! strcmp(pdb.resAtoms[j].names[k], name))
					{
						for (l=0; l<nbTypeRad; l++)
						{
							if (! strcmp(pdb.resAtoms[j].types[k], radia[l].type))
							{
								rad[i] = radia[l].radius;
								break;
							}
						}
						break;
					}
				}
				break;
			}
		}

		// If, in the database, there is no residue, or no atom with the name specified in the PDB file,
		// or if the residue contains atom specific of termini, then rad[i] = 0.
		if (rad[i] < dEpsilon)
		{
			if (termres >= 0)	// Testing whether the residue contains termini atoms
			{
				for (k=0; k<pdb.resAtoms[termres].nbAtoms; k++)
				{
					if (! strcmp(pdb.resAtoms[termres].names[k], name))
					{
						for (l=0; l<nbTypeRad; l++)
						{
							if (! strcmp(pdb.resAtoms[termres].types[k], radia[l].type))
							{
								rad[i] = radia[l].radius;
								break;
							}
						}
						break;
					}
				}
			}
		}

		// This time, if there is no termini atoms, the database does not contain the required information
		if (rad[i] < dEpsilon)
		{
			strncpy(resi, &(*ENSData).atomHeaders[i][22], 5);
			deepTrim(resi);

			printf("ERROR: Impossible to attribute a radius to atom %s of residue %s %s!\n", name, resn, resi);
			_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
			return access;
		}
	}


	// Calculate the lowest zbound for each actom and sort them from low to high on zlb
	if (radtype == 'S') radFact = .5;		// Half the sigma value
	for (i=0; i<dim; i++)
	{
		rad[i] *= radFact;
		zlb[i] = (*ENSData).atomCoords[0][i].z - rad[i];
	}

	sortag(zlb, tag, dim);		// Sort zlb and impose the same order to tag

	for (i=0; i<dim; i++)
	{
		j = tag[i];
		tag1[j] = i;
		rad1[i] = rad[j];

		xr1[i] = (*ENSData).atomCoords[0][j].x;
		yr1[i] = (*ENSData).atomCoords[0][j].y;
		zr1[i] = (*ENSData).atomCoords[0][j].z;

		if (rad1[i] < rmin) rmin = rad1[i];
	}

	// The radius of an atom sphere = atom radius + probe radius
	zor = zlb[0] - rh2o;
	for (i=0; i<dim; i++)
	{
		rad1[i] += rh2o;
		zr1[i]  -= zor;			// The atom that had the lowest Z coordinate has now a Z coordinate of rh2o
		zub[i]   = zr1[i] + rad1[i];
		zlb[i]   = zr1[i] - rad1[i];
		if (zub[i] > zubmax) zubmax = zub[i];
		area[i] = 0.;
	}

	// Z resolution determined
	hzres = (rmin + rh2o) * accuracy;
	access.zres = 2. * hzres;
	cutoff = access.zres / 100.;
	iano = 0;		// Index of atoms in tables
	j1 = 0;
	iflag = 0;
	znext = access.zres;
	access.zsecnb = zubmax / access.zres;
	access.icnt1  = 0;
	access.icnt2  = 0;

	// Section atom spheres perpendicular to the z axis
	for (i=0; i<access.zsecnb; i++)
	{
		itab = 0;
		zgrid = znext;
		znext += access.zres;

		for (j=iano; j<dim; j++)
		{
			// The upper and lower z bounds of an atom are used to determine
			// whether this z-section cuts the atom sphere
			if (zub[j] > zgrid)
			{
				if (zlb[j] >= zgrid) break;
				
				kn[j] = -1;
				for (k=0; k<ict; k++) arci[k][j] = 0.;

				// Count and record spheres cut by section
				inz[itab] = j;
				++itab;

				// Find radius of circle locus
				a = rad1[j] * rad1[j] - pow(zgrid-zr1[j], 2);
				rsec2[j] = a;
				rsec[j] = sqrt(a);
			}
			
			if (iflag == 1) continue;

			// Find first sphere cut by section I, skip atoms preceeding it in list for next section
			if (zub[j] > znext)
			{
				iflag = 1;
			}
			else
			{
				++j1;
			}
		}
		iano = j1;
		iflag = 0;

		// Zero, one and more circles on section require different processing
		if (itab < 1) continue;
		if (itab != 1)
		{
			//j3 = itab - 1;

			// Find intersections of circles in section calarc called to find
			// initial and final angles of intersection of circles. If arci and
			// arcf arrays are filled, reduce is called.
			// The current first index of these arrays for the atom indicated by 
			// the second index is stored in the array kn. If kn for any atom is
			// 10000, then no area remains accessible for that atom on this section
			// or the atom is not of interest (knset=1)
			for (j=0; j<itab-1; j++)
			{
				l = inz[j];			// Index of atom 1 that is cut by the z-section
				for (k=j+1; k<itab; k++)
				{
					m = inz[k];		// Index of atom 2 that is cut by the z-section
					a = rsec[l] + rsec[m];

					dx = xr1[m] - xr1[l];
					if (fabs(dx)>=a) continue;

					dy = yr1[m] - yr1[l];
					if (fabs(dy)>=a) continue;

					d2 = dx * dx + dy * dy;
					d = sqrt(d2);
					if (d>=a) continue;

					if (kn[l] >= ict1-1 && kn[l] < ict)
					{
						reduce(&l, 0, ict, arci, arcf, tag, tag1, zr1, rad1, area, inz, kn, &access.icnt1, &access.icnt2, ict1, access.zres, rh2o, zgrid , hzres, cutoff);
					}
					++kn[l];

					if (kn[m] >= ict1-1 && kn[m] < ict)
					{
						reduce(&m, 0, ict, arci, arcf, tag, tag1, zr1, rad1, area, inz, kn, &access.icnt1, &access.icnt2, ict1, access.zres, rh2o, zgrid , hzres, cutoff);
					}
					++kn[m];

					// Do the circles intersect, or is one completely inside the other ?
					b = rsec[m] - rsec[l];

					if (d > fabs(b))
					{
						q = rsec2[m] - rsec2[l];
						d *= 2.;
						if (kn[m] <= ict1) calarc(m,  1, d2, q, d, rsec, arci, arcf, kn, dy, dx);
						if (kn[l] <= ict1) calarc(l, -1, d2, q, d, rsec, arci, arcf, kn, dy, dx);
					}
					else if (b > 0.)
					{
						kn[l] = 9999;
						--kn[m];
					}
					else
					{
						kn[m] = 9999;
						--kn[l];
					}
					
				}
			}
		}

		// Find the accessible contact surface area for all spheres intersecting this section
		reduce(&idum, itab, ict, arci, arcf, tag, tag1, zr1, rad1, area, inz, kn, &access.icnt1, &access.icnt2, ict1, access.zres, rh2o, zgrid , hzres, cutoff);
	}

	// Output operation parameters
	if (mode == 'A')
	{
		for (i=0; i<dim; i++)
		{
			j = tag1[i];
			access.as[i] = area[j];
		}
	}
	else
	{
		for (i=0; i<dim; i++)
		{
			j = tag1[i];
			access.as[i] = area[j] * pow((rad1[j]-rh2o)/rad1[j], 2);
		}
	}

	_cleanh(ict, tag, tag1, kn, inz, zlb, zub, rad, rad1, xr1, yr1, zr1, rsec, rsec2, area, arci, arcf, pdb, radia);
	access.FAILED = 0;
	return access;
}


#endif