// 01/25/2007, fix the probelm of missed ring atoms for ring current shifts calculation
#include "PDB.h"

PDB::PDB()
{
	r1 = 9999; rN = -9999;
}


PDB::PDB(const string &fileName)
{
	PDBfileName = fileName;

	r1 = 9999; rN = -9999;

	loadPDB(fileName);
}


string PDB::getField(const string &str, int index)
{
	switch(index)
	{
		case 1: return str.substr(0,6); //"ATOM"
		case 2: return str.substr(6,5); //Atom number
		case 3: return str.substr(11,5); //Atom name
		case 4: return str.substr(17,4); //Residue name
		case 5: return str.substr(21,1); //Chain ID
		case 6: return str.substr(22,4); //Reside seq
		case 7: return str.substr(30,8); //X
		case 8: return str.substr(38,8); //Y
		case 9: return str.substr(46,8); //Z
		case 10: return str.substr(54,6); //Occupancy
		case 11: return str.substr(60,6); // B-factor
	}

	return NULL;
}


void PDB::loadPDB_Entry(const string &str, PDB_Entry &entry)
{
	entry.atomNum = atoi( getField(str,2).c_str() );
	
	string atomName = simplifyWhiteSpace( getField(str,3).c_str() );
	entry.atomName = atomName;
	if( atomName == "H" ) entry.atomName="HN";
	if( atomName[0] >= '1' && atomName[0] <= '3' ) entry.atomName = atomName.substr(1,4) + atomName[0];
	
	// remove the charge character '+' 
	entry.resName =  simplifyWhiteSpace( getField(str,4).substr(0,3).c_str() );
	entry.chainName = getField(str,5);
	entry.resNum = atoi( getField(str,6).c_str() );
	entry.X = atof( getField(str,7).c_str() );
	entry.Y = atof( getField(str,8).c_str() );
	entry.Z = atof( getField(str,9).c_str() );
	entry.B_Factor = atof( getField(str,11).c_str() );;

	entry.Coord[0] = entry.X;
	entry.Coord[1] = entry.Y;
	entry.Coord[2] = entry.Z;
}


void PDB::loadPDB(const string &fileName)
{
	PDBfileName = fileName;

	ifstream file(fileName.c_str());
	if(! file.is_open() )
	{
		cerr << "\tCan't open file " << fileName << " for reading" << endl;
		exit(0);
	}

	residList.clear();
	residListOne.clear();
	ATOMS.clear();
	Conformers.clear();
	
	acceptorList.clear();
	donorList.clear();
	HBDistList.clear();
	

	Mol conformer;
	string str;
	int conformer_count=1;

	while (! file.eof() ) 
	{
		getline(file, str);
		if(str.empty() || str.size() == 0) continue;

		int pos = str.find("ATOM"); //starting of one molecule
		int pos2 = str.find("TER"), pos3 = str.find("END"); //ending of one molecule
		if(pos != 0 && pos2 != 0 && pos3 != 0) continue;
		
		if( pos2 == 0 || pos3 == 0) 
		{
			Conformers[conformer_count] = conformer;
			conformer_count++;
			conformer.clear();
			break; // ******* only read one confomer/chain
			//continue;
		}

		PDB_Entry temp;
		loadPDB_Entry(str, temp);

		conformer[temp.atomNum] = temp; // table indexed by the atom number
		ATOMS[conformer_count][temp.resNum][temp.atomName] = temp; 
		// table indexed by the residue number and atom name
		// fast for obtaining a specific atom with given residue number and atom name

		residList[temp.resNum] = temp.resName; //resdiue list with three-letter-aa name
		residListOne[temp.resNum] = getOneAAName(temp.resName); //resdiue list with one-letter-aa name

		if( temp.resNum < r1 ) r1 = temp.resNum;
		if( temp.resNum > rN ) rN = temp.resNum;

		//sprintf(buf, "%5d %-4s %3s %1s %3d    %8.3f%8.3f%8.3f %6.2f",temp.atomNum,(const char*)(temp.atomName),(const char*)(temp.resName),\
		//	(const char*)(temp.chainName),temp.resNum,temp.X,temp.Y,temp.Z,temp.B_Factor);

		//cout << temp.resNum << endl ;
	}
	file.close();

	if(conformer.size() != 0) Conformers[conformer_count] = conformer;

	//correction of the name cys/CYS accoding to the status of the disulfide bond
	map<int, string> ::iterator it;
	for(it=residListOne.begin(); it!=residListOne.end(); it++)
		if( it->second == "C" )
			if( isSSBonded(1, it->first) ) {
				residListOne[it->first] = "c";
				residList[it->first] = "cys";
			}

}


bool PDB::isSSBonded(int conformerID, int resNum)
{
	if( residListOne[resNum] != "C" && residListOne[resNum] != "c" ) return false;

	PDB_Entry CYS_SG = getEntry(1,resNum, "SG");

	map<int, string>::iterator it;
	for ( it = residListOne.begin(); it != residListOne.end(); it++ )
	{
		if( (it->second == "C" || it->second == "c") && abs(it->first - resNum) >= 4 )
			if( getDist(CYS_SG.Coord, ATOMS[1][it->first]["SG"].Coord ) <= 2.5 )
				return true;
	}

	return false;
}



PDB_Entry PDB::getEntry(int conformerID, int rNum, const string &aName)
{
	return ATOMS[conformerID][rNum][aName];
}



PDB_Entry PDB::getEntry(int conformerID, int aNum)
{
	return Conformers[conformerID][aNum]; 
}


float PDB::getBondAngle(Vec3 A, Vec3 B, Vec3 C)
{
	Vec3 v,v1,v2;

	Vec3Copy(v,B); Vec3Copy(v1,A); Vec3Copy(v2,C);
	Vec3Sub(v1,v); // vector BA
	Vec3Sub(v2,v); // vector BC

	Vec3Norm(v1); // unit vector BA
	Vec3Norm(v2); // unit vector BC

	float c = Vec3Scalar(v1, v2); // scalar products
	Vec3Cross(v1, v2); // cross product
	float s = Vec3Abs(v1);
  
	return atan2f(s, c)*RAD;
}


float PDB::getBondAngle(PDB_Entry a, PDB_Entry b, PDB_Entry c)
{  
	return getBondAngle(a.Coord, b.Coord, c.Coord);
}


float PDB::getDihedralAngle(PDB_Entry a, PDB_Entry b, PDB_Entry c, PDB_Entry d)
{
	double cb[3], n1[3], n2[3];
	float co;
	double TEMP, TEMP1, TEMP2, TEMP3, TEMP4, TEMP5;

	if( (a.atomName).empty() || (b.atomName).empty() || (c.atomName).empty() || (d.atomName).empty() ) 
		return MAXNUM;

	cb[0] = c.X - b.X; cb[1] = c.Y - b.Y; cb[2] = c.Z - b.Z; // vector b->c
	n1[0] = (b.Y - a.Y) * cb[2] + (a.Z - b.Z) * cb[1]; //normal vector of plane a-b-c
	n1[1] = (b.Z - a.Z) * cb[0] + (a.X - b.X) * cb[2];
	n1[2] = (b.X - a.X) * cb[1] + (a.Y - b.Y) * cb[0];
	n2[0] = cb[1] * (d.Z - c.Z) + cb[2] * (c.Y - d.Y); //normal vector of plane b-c-d
	n2[1] = cb[2] * (d.X - c.X) + cb[0] * (c.Z - d.Z);
	n2[2] = cb[0] * (d.Y - c.Y) + cb[1] * (c.X - d.X);
	TEMP  = n1[0]; TEMP1 = n1[1]; TEMP2 = n1[2];
	TEMP3 = n2[0]; TEMP4 = n2[1]; TEMP5 = n2[2];
	co = (float)((n1[0] * n2[0] + n1[1] * n2[1] + n1[2] * n2[2]) /
	       sqrt( (TEMP*TEMP + TEMP1*TEMP1 + TEMP2*TEMP2)*(TEMP3*TEMP3 + TEMP4*TEMP4 + TEMP5*TEMP5) ) );
	// cos(a-b-c-d) = |n1.n2|/|n1|.|n2|       
	
	return RAD*(arccos_(co)*
		sgn( (float)( (n1[1]*n2[2] - n1[2]*n2[1])*cb[0] + (n1[2]*n2[0] - n1[0]*n2[2])*cb[1] + (n1[0]*n2[1] - n1[1]*n2[0])*cb[2]) ) );
}


float PDB::getPhi(int conformerID, int resNum)
{
	return getDihedralAngle(getEntry(conformerID,resNum-1,"C"), \
		getEntry(conformerID,resNum,"N"), \
		getEntry(conformerID,resNum,"CA"), \
		getEntry(conformerID,resNum,"C") );
}


float PDB::getPsi(int conformerID, int resNum)
{
	return getDihedralAngle(getEntry(conformerID,resNum,"N"), \
		getEntry(conformerID,resNum,"CA"), \
		getEntry(conformerID,resNum,"C"), \
		getEntry(conformerID,resNum+1,"N") );
}


float PDB::getOmega(int conformerID, int resNum)
{
	return getDihedralAngle(getEntry(conformerID,resNum-1,"CA"), \
		getEntry(conformerID,resNum-1,"C"), \
		getEntry(conformerID,resNum,"N"), \
		getEntry(conformerID,resNum,"CA"));
}


float PDB::getChi1(int conformerID, int resNum)
{
	string rName = residList[resNum];

	PDB_Entry R;
	if( rName == "GLY" || rName == "ALA" ) 
		return MAXNUM;
	else if( rName == "VAL" )
		R = getEntry(conformerID,resNum,"CG2");
	else if( rName == "ILE" ) 
		R = getEntry(conformerID,resNum,"CG1");
	else if( rName == "THR" ) 
		R = getEntry(conformerID,resNum,"OG1");
	else if( rName == "SER" ) 
		R = getEntry(conformerID,resNum,"OG" );
	else if( rName == "CYS" || rName == "cys") 
		R = getEntry(conformerID,resNum,"SG" );
	else
		R = getEntry(conformerID,resNum,"CG" );

	return getDihedralAngle(getEntry(conformerID,resNum,"N"), \
							getEntry(conformerID,resNum,"CA"), \
							getEntry(conformerID,resNum,"CB"), R);
}


float PDB::getChi2(int conformerID, int resNum)
{
	string rName = residList[resNum];
	if( rName.compare("GLY") == 0 || rName.compare("ALA") == 0) 
		return MAXNUM;

	PDB_Entry R3, R4;

	if( rName == "VAL" ){ //chi21
		R3 = getEntry(conformerID,resNum,"CG1");
		R4 = getEntry(conformerID,resNum,"HG11");
	}
	else if( rName == "THR" ) {
		R3 = getEntry(conformerID,resNum,"CG2");
		R4 = getEntry(conformerID,resNum,"HG21");
	}
	else if( rName == "SER" ) {
		R3 = getEntry(conformerID,resNum,"OG");
		R4 = getEntry(conformerID,resNum,"HG");
	}
	else if( rName == "ILE" ) { //chi21
		R3 = getEntry(conformerID,resNum,"CG1" );
		R4 = getEntry(conformerID,resNum,"CD1" );
	}
	else if( rName == "CYS" || rName == "cys") {
		R3 = getEntry(conformerID,resNum,"SG" );
		R4 = getEntry(conformerID,resNum,"HG" );
	}
	else
		R3 = getEntry(conformerID,resNum,"CG" );

	if( rName == "ASN" || rName == "ASP" ) R4 = getEntry(conformerID,resNum,"OD1" );
	else if( rName == "HIS" ) R4 = getEntry(conformerID,resNum,"ND1" );
	else if( rName == "MET" ) R4 = getEntry(conformerID,resNum,"SD" );
	else if( rName == "LEU" || rName == "PHE" || rName == "TYR" || rName == "TRP" ) 
		R4 = getEntry(conformerID,resNum,"CD1" );
	else R4 = getEntry(conformerID,resNum,"CD" );

	return getDihedralAngle(getEntry(conformerID,resNum,"CA"),getEntry(conformerID,resNum,"CB"),R3, R4);
}


string PDB::getThreeAAName(char a)
{
	switch(a)
	{
		case 'A': return "ALA";
		case 'C': return "CYS";
		case 'c': return "cys";
		case 'D': return "ASP";
		case 'E': return "GLU";
		case 'F': return "PHE";
		case 'G': return "GLY";
		case 'H': return "HIS";
		case 'I': return "ILE";
		case 'K': return "LYS";
		case 'L': return "LEU";
		case 'M': return "MET";
		case 'N': return "ASN";
		case 'P': return "PRO";
		case 'Q': return "GLN";
		case 'R': return "ARG";
		case 'S': return "SER";
		case 'T': return "THR";
		case 'V': return "VAL";
		case 'W': return "TRP";
		case 'Y': return "TYR";
	}

	return "???";
}



string PDB::getOneAAName(const string& a)
{
	if( a == "ALA") return "A";
	else if( a == "CYS") return "C";
	else if( a == "cys") return "c";
	else if( a == "ASP") return "D";
	else if( a == "GLU") return "E";
	else if( a == "PHE") return "F";
	else if( a == "GLY") return "G";
	else if( a == "HIS" || a == "HIS+") return "H";
	else if( a == "ILE") return "I";
	else if( a == "LYS" || a == "LYS+") return "K";
	else if( a == "LEU") return "L";
	else if( a == "MET") return "M";
	else if( a == "ASN") return "N";
	else if( a == "PRO") return "P";
	else if( a == "GLN") return "Q";
	else if( a == "ARG" || a == "ARG+") return "R";
	else if( a == "SER") return "S";
	else if( a == "THR") return "T";
	else if( a == "VAL") return "V";
	else if( a == "TRP") return "W";
	else if( a == "TYR") return "Y";

	return "?";
}


//collet the ring information required for calculating Haigh-Mallion ring shifts
//Haigh, C.W. and Mallion, R.B. (1979) Progr. NMR Spectrosc., 13, 303-344
//Parameters are taken from Case, D.A. (1995) J. Biomol.NMR, 6, 341-346.
void PDB::initOrbitalShift(){

	Mol conf = Conformers[1];

	string PF6Atoms[6] = {"CG", "CD1", "CE1", "CZ", "CE2", "CD2" };
	string W6Atoms[6] = {"CD2", "CE2", "CZ2", "CH2", "CZ3", "CE3" };
	string H5Atoms[5] = {"CG", "ND1", "CE1", "NE2", "CD2" };
	string W5Atoms[5] = {"CG", "CD1", "NE1", "CE2", "CD2" };

	map<int, string>::iterator it;
	Mol::iterator itA;

	int i;
	RingNo = 0;

	for ( it = residList.begin(); it != residList.end(); it++ )
	{
		int resID = it->first; //.resNum;
		string resName = it->second;
		
		if( resName == "PHE" || resName == "TYR" || resName == "TRP") // 6-member ring
		{
			Rings[RingNo].resID = resID; Rings[RingNo].resName = resName; Rings[RingNo].atomNo = 6;
			if( resName == "PHE" ) Rings[RingNo].ringFact = 1.46;
			if( resName == "TYR" ) Rings[RingNo].ringFact = 1.24;
			if( resName == "TRP" ) Rings[RingNo].ringFact = 1.24;

			bool NO_MISSED_RING_ATOMS = true;

			for( i = 0; i < 6; i++ )
			{
				PDB_Entry atom;
				if( resName == "PHE" || resName == "TYR")
				{
					atom = getEntry(1,resID,PF6Atoms[i]);
					if( atom.atomName == PF6Atoms[i] )
						Vec3Copy(Rings[RingNo].coordA[i], atom.Coord);
					else NO_MISSED_RING_ATOMS = false; // missed ring atom in coordinates file
				}
				else if ( resName == "TRP")
				{
					atom = getEntry(1,resID,W6Atoms[i]);
					if( atom.atomName == W6Atoms[i] )
						Vec3Copy(Rings[RingNo].coordA[i], atom.Coord);
					else NO_MISSED_RING_ATOMS = false; // missed ring atom in coordinates file
				}
			}

			if(NO_MISSED_RING_ATOMS) RingNo++;
		}
		
		if( resName == "HIS" || resName == "TRP")
		{ // 5-member ring
			Rings[RingNo].resID = resID; Rings[RingNo].resName = resName; Rings[RingNo].atomNo = 5;
			if( resName == "HIS" ) Rings[RingNo].ringFact = 1.35;
			if( resName == "TRP" ) Rings[RingNo].ringFact = 1.32;

			bool NO_MISSED_RING_ATOMS = true;

			for( i = 0; i < 5; i++ )
			{
				PDB_Entry atom;
				if( resName == "HIS" )
				{
					atom = getEntry(1,resID,H5Atoms[i]);
					if( atom.atomName == H5Atoms[i] )
						Vec3Copy(Rings[RingNo].coordA[i], atom.Coord);
					else NO_MISSED_RING_ATOMS = false; // missed ring atom in coordinates file
				}
				else if ( resName == "TRP" )
				{
					atom = getEntry(1,resID,W5Atoms[i]);
					if( atom.atomName == W5Atoms[i] )
						Vec3Copy(Rings[RingNo].coordA[i], atom.Coord);
					else NO_MISSED_RING_ATOMS = false; // missed ring atom in coordinates file
				}
			}

			if(NO_MISSED_RING_ATOMS) RingNo++;
		}
	}

	for(i=0; i < RingNo; i++)
		calcPlane(&Rings[i]);
}


//Calculate the parameters for a ring
void PDB::calcPlane(RingData *ringP)
{
	Vec3 v1, v2;
	int atomI;

	Vec3Zero(ringP->center);
	for(atomI = 0; atomI < ringP->atomNo; atomI++)
		Vec3Add(ringP->center, ringP->coordA[atomI]);
	Vec3Scale(ringP->center, 1.0f/ringP->atomNo);

	Vec3Zero(ringP->norm);
	ringP->rad = 0.0f;
	Vec3Copy(v1, ringP->coordA[ringP->atomNo - 1]);
	Vec3Sub(v1, ringP->center);
	for(atomI = 0; atomI < ringP->atomNo; atomI++) 
	{
		Vec3Copy(v2, ringP->coordA[atomI]);
		Vec3Sub(v2, ringP->center);
		ringP->rad += Vec3Abs(v2);
		Vec3Cross(v1, v2);
		Vec3Add(ringP->norm, v1);
		Vec3Copy(v1, v2);
	}
	Vec3Norm(ringP->norm);
	ringP->rad /= ringP->atomNo;

	/* transformation matrix ring plane -> x-y plane */
	ringP->transM[0][2] = ringP->norm[0];
	ringP->transM[1][2] = ringP->norm[1];
	ringP->transM[2][2] = ringP->norm[2];

	 /* vector orthgonal to norm */
	if (ringP->norm[0] > 0.5f || ringP->norm[0] < -0.5f) 
	{
		v1[0] = - ringP->norm[1];
		v1[1] = ringP->norm[0];
		v1[2] = 0.0f;
	} else 
	{
		v1[0] = 0.0f;
		v1[1] = - ringP->norm[2];
		v1[2] = ringP->norm[1];
	}
	Vec3Norm(v1);
	ringP->transM[0][1] = v1[0];
	ringP->transM[1][1] = v1[1];
	ringP->transM[2][1] = v1[2];

	Vec3Cross(v1, ringP->norm);
	ringP->transM[0][0] = v1[0];
	ringP->transM[1][0] = v1[1];
	ringP->transM[2][0] = v1[2];
}



float PDB::getOrbitalShift(int conformerID, int resNum, const string &aName)
{
	Vec3 v;

	Vec3Copy( v, getEntry(conformerID,resNum,aName).Coord );

	//cout << RingNo << " " << resNum << " " << aName << endl;

	float ringShiftSum = 0.0f;
	//loop over all rings in protein
	for(int i=0; i < RingNo; i++) 
	{
		int atomNo, atomI;
		Vec3 vt, v1, v2;
		float area, r1, r2, g, b;

		//exclude the effects from aromatic rings themselves (except for HN atom)
		if( Rings[i].resID == resNum && aName != "HN" ) continue;
		
		atomNo = Rings[i].atomNo;
		g = 0.0f;

		//loop over all 6(or 5) atoms in a ring
		for (atomI = 0; atomI < atomNo; atomI++) 
		{
			Vec3Copy(v1, Rings[i].coordA[atomI]);
			Vec3Copy(v2, Rings[i].coordA[(atomI + 1) % atomNo]);

			r1 = Vec3DiffAbs(v1, v);
			r2 = Vec3DiffAbs(v2, v);

			Vec3Copy(vt, v);
			Vec3Sub(vt, Rings[i].center);
			Vec3Sub(v1, Rings[i].center);
			Vec3Sub(v2, Rings[i].center);

			Mat3VecMult(vt, Rings[i].transM);
			Mat3VecMult(v1, Rings[i].transM);
			Mat3VecMult(v2, Rings[i].transM);
			Vec3Sub(v1, vt);
			Vec3Sub(v2, vt);

			area = v1[0] * v2[1] - v1[1] * v2[0];
			g += area * (1.0f / (r1 * r1 * r1) + 1.0f / (r2 * r2 * r2));
			//cout << atomI << " " << g << " ";
		}

		g *= 0.5f;
		b = (float) 5.4548e-6;

		ringShiftSum += Rings[i].ringFact * b * g;

		//if( Rings[i].ringFact * b * g *((float)-1.0e6) > 0.6)
		//	cout << ringShiftSum << " " << g << " " << Rings[i].resID  << resNum << " " << aName << endl;
	}

	//cout << ringShiftSum*((float)-1.0e6) << " " << resNum << " " << aName << endl;

	return ringShiftSum*((float)-1.0e6);
}



float PDB::getDist(PDB_Entry A, PDB_Entry B)
{
    return getDist(A.Coord, B.Coord);
}



float PDB::getDist(Vec3 A, Vec3 B)
{
	Vec3 v;

    Vec3Copy(v, B);
    Vec3Sub(v, A);

    return Vec3Abs(v);
}



void PDB::initHBond(float DIST, float ANGLE)
{
	Mol conf = Conformers[1];
	Mol::iterator it;

	acceptorList.clear();
	donorList.clear();
	HBDistList.clear();

	//retreive the donor and acceptor list from coordinates
	for ( it = conf.begin(); it != conf.end(); it++ )
	{
		PDB_Entry a = isAcceptor(it->second);

		if( !(a.atomName.empty()) ) acceptorList[it->first] = a.atomNum;
		else {
			PDB_Entry temp = isDonor(it->second);
			if(temp.atomName.empty()) continue;
			donorList[it->first] = temp.atomNum;
		}
	}


	map<int, int>::iterator itA, itD;
	for(itA = acceptorList.begin(); itA != acceptorList.end(); itA++ ) //loop over acceptor list
	{
		PDB_Entry A = conf[ itA->first ], A_Heavy = conf[ itA->second ];;
		
		//search for donors
		for( itD = donorList.begin(); itD != donorList.end(); itD++ )
		{
			PDB_Entry D = conf[ itD->first ], D_Heavy = conf[ itD->second ];

			if( abs( A.resNum - D.resNum) < 2) continue; //minimal 2 residues apart

			float D_ON = getDist(A,D_Heavy);
			float D_CH = getDist(A_Heavy,D);
			float D_OH = getDist(A,D);
			float D_CN = getDist(A_Heavy,D_Heavy);
			
			if( D_OH > D_CN ) continue;

			float HBond_E = 332.0*0.42*0.20;
			HBond_E *= ( 1.0/(D_ON*D_ON*D_ON)+1.0/(D_CH*D_CH*D_CH)-1.0/(D_OH*D_OH*D_OH)-1.0/(D_CN*D_CN*D_CN) );
			//Kabsch, W. and Sander, C. (1983) Biopolymer, 22, 2577-2637.

			if( HBond_E < -0.5) //is a HB bond
			{
				if( HBDistList[A.resNum][A.atomName] == 0 || HBDistList[A.resNum][A.atomName] > D_OH ) 
					HBDistList[A.resNum][A.atomName] = D_OH;
				// keep the strongest Hbond

				if( HBDistList[D.resNum][D.atomName] == 0 || HBDistList[D.resNum][D.atomName] > D_OH) 
					HBDistList[D.resNum][D.atomName] = D_OH;
				// keep the strongest Hbond
			}
		}
	}

}


// must run initHBond() first
float PDB::getHBondDist(PDB_Entry D)
{
	return HBDistList[D.resNum][D.atomName]; 
}


// must run initHBond() first
float PDB::getHBondDist(int resNum, string atomName)
{
	return HBDistList[resNum][atomName]; 
}



PDB_Entry PDB::isAcceptor(PDB_Entry A)
{
	if(A.atomName == "O" || A.atomName == "OT1"  || A.atomName == "OT2" ) return ATOMS[1][A.resNum]["C"];

	if(A.resName == "ASN" && A.atomName == "OD1") return ATOMS[1][A.resNum]["CG"];
	if(A.resName == "ASP" && (A.atomName == "OD1" || A.atomName == "OD2") ) return ATOMS[1][A.resNum]["CG"];
	if((A.resName == "CYS"||A.resName == "cys") && A.atomName == "SG") return ATOMS[1][A.resNum]["CB"];
	if(A.resName == "GLN" && A.atomName == "OE1") return ATOMS[1][A.resNum]["CD"];
	if(A.resName == "GLU" && (A.atomName == "OE1" || A.atomName == "OE2") ) return ATOMS[1][A.resNum]["CD"];
	if(A.resName == "HIS" && (A.atomName == "ND1" || A.atomName == "NE2") ) return ATOMS[1][A.resNum]["CE1"]; // could be CG or CD2
	if(A.resName == "MET" && A.atomName == "SD") return ATOMS[1][A.resNum]["CG"];
	if(A.resName == "SER" && A.atomName == "OG") return ATOMS[1][A.resNum]["CB"];
	if(A.resName == "THR" && A.atomName == "OG1") return ATOMS[1][A.resNum]["CB"];
	if(A.resName == "TYR" && A.atomName == "OH") return ATOMS[1][A.resNum]["CZ"];

	return EMPTY;
}


PDB_Entry PDB::isDonor(PDB_Entry D)
{
	if(D.atomName == "HN") return ATOMS[1][D.resNum]["N"];

	if(D.atomName.find("HA") != string::npos ) return ATOMS[1][D.resNum]["CA"];

	if(D.resName == "ARG") {
		if(D.atomName == "HE") return ATOMS[1][D.resNum]["NE"];
		if(D.atomName == "HH11" || D.atomName == "HH12") return ATOMS[1][D.resNum]["NH1"];
		if(D.atomName == "HH21" || D.atomName == "HH22") return ATOMS[1][D.resNum]["NH2"];
	}
	if(D.resName == "ASP" && (D.atomName == "HD21" || D.atomName == "HD22") ) return ATOMS[1][D.resNum]["ND2"];
	if( (D.resName == "CYS"||D.resName == "cys") && D.atomName == "HG" ) return ATOMS[1][D.resNum]["SG"];
	if(D.resName == "GLU" && (D.atomName == "HE21" || D.atomName == "HE22") ) return ATOMS[1][D.resNum]["NE2"];
	if(D.resName == "HIS") {
		if(D.atomName == "HD1") return ATOMS[1][D.resNum]["ND1"];
		if(D.atomName == "HE2") return ATOMS[1][D.resNum]["NE2"];
	}
	if(D.resName == "LYS" && (D.atomName == "HZ1" || D.atomName == "HZ2" || D.atomName == "HZ3") ) return ATOMS[1][D.resNum]["NZ"];
	if(D.resName == "SER" && D.atomName == "HG" ) return ATOMS[1][D.resNum]["OG"];
	if(D.resName == "THR" && D.atomName == "HG1" ) return ATOMS[1][D.resNum]["OG1"];
	if(D.resName == "TRP" && D.atomName == "HE1" ) return ATOMS[1][D.resNum]["NE1"];
	if(D.resName == "TYR" && D.atomName == "HH" ) return ATOMS[1][D.resNum]["OH"];

	return EMPTY;
}



long PDB::sgn(float x)
{
	return ((x >= 0.0) - (x < 0.0));
}


float PDB::arccos_(float x)
{
	if (x > 1.0) x = 1.0;
	else if (x < -1.0) x = -1.0;

	if (fabs(x) >= 0.5) {
		if (x > 0.0) 
			return atan(sqrt(1.0 - x * x) / x);
		else 
			return (PI + atan(sqrt(1.0 - x * x) / x));
	} 
	else
		return (0.5 * PI - atan(x / sqrt(1.0 - x * x)));
}


void PDB::Vec3Zero(Vec3 v)
{
	for (int i = 0; i < 3; i++)
		v[i] = 0.0f;
}


void PDB::Vec3Copy(Vec3 v1, Vec3 v2)
{
	for (int i = 0; i < 3; i++)
		v1[i] = v2[i];
}


float PDB::Vec3Abs(Vec3 v)
{
	return sqrtf(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}


float PDB::Vec3DiffAbs(Vec3 v1, Vec3 v2)
{
	return sqrtf((v1[0] - v2[0]) * (v1[0] - v2[0]) + 
		(v1[1] - v2[1]) * (v1[1] - v2[1]) +
		(v1[2] - v2[2]) * (v1[2] - v2[2]));
}


void PDB::Vec3Norm(Vec3 v)
{
	float a = Vec3Abs(v);
	for(int i = 0; i < 3; i++)
		v[i] /= a;
}


void PDB::Vec3Scale(Vec3 v, float s)
{
	for(int i = 0; i < 3; i++)
		v[i] *= s;
}


void PDB::Vec3Add(Vec3 v1, Vec3 v2)
{
	for(int i = 0; i < 3; i++)
		v1[i] += v2[i];
}


void PDB::Vec3Sub(Vec3 v1, Vec3 v2)
{
	for(int i = 0; i < 3; i++)
		v1[i] -= v2[i];
}


float PDB::Vec3Scalar(Vec3 v1, Vec3 v2)
{
	return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}


void PDB::Vec3Cross(Vec3 v1, Vec3 v2)
{
	Vec3 vRes;

	vRes[0] = v1[1] * v2[2] - v1[2] * v2[1];
	vRes[1] = v1[2] * v2[0] - v1[0] * v2[2];
	vRes[2] = v1[0] * v2[1] - v1[1] * v2[0];
	Vec3Copy(v1, vRes);
}


void PDB::Mat3VecMult(Vec3 v, Mat3 m)
{
	Vec3 vRes;
 
	vRes[0] = v[0] * m[0][0] + v[1] * m[1][0] + v[2] * m[2][0];
	vRes[1] = v[0] * m[0][1] + v[1] * m[1][1] + v[2] * m[2][1];
	vRes[2] = v[0] * m[0][2] + v[1] * m[1][2] + v[2] * m[2][2];
 
	for(int i = 0; i < 3; i++)
		v[i] = vRes[i];
}


string PDB::simplifyWhiteSpace(const string &str)
{
	if ( str.empty() )				// nothing to do
		return str;
    
	string result;
	result.resize(str.length()); //.setLength( length() );
	const char *from = str.c_str();
	const char *fromend = from+str.length();
	int outc=0;
	char *to	= &(result[0]);

	while( true ) 
	{
		while( from!=fromend && isSpace(*from) )
			from++;
		while( from!=fromend && !isSpace(*from) )
			to[outc++] = *from++;
		if( from!=fromend )
			to[outc++] = ' ';
		else
			break;
	}
    
	if ( outc > 0 && to[outc-1] == ' ' )
		outc--;
	result.resize( outc );
	return result;
}


bool PDB::isSpace( const char &c )
{
	if( (c >= 9 && c <= 13) || c == ' ' ) return true;
	return false;
}




