// 01/25/2007, fix the probelm of missed ring atoms for ring current shifts calculation #include "PDB.h" PDB::PDB() { Format = "%4s %5d %-4s %3s %1s %3d %8.3f%8.3f%8.3f%6.2f%6.2f"; VarsNumber = 11; r1 = 9999; rN = -9999; } PDB::PDB(const string &fileName) { PDBfileName = fileName; Format = "%4s %5d %-4s %3s %1s %3d %8.3f%8.3f%8.3f%6.2f%6.2f"; VarsNumber = 11; 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() ); entry.atomName = simplifyWhiteSpace( getField(str,3).c_str() ); if( entry.atomName == "H" ) entry.atomName="HN"; // 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(); 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 ::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::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); Vec3Sub(v2, v); Vec3Norm(v1); Vec3Norm(v2); float c = Vec3Scalar(v1, v2); Vec3Cross(v1, v2); 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; n1[0] = (b.Y - a.Y) * cb[2] + (a.Z - b.Z) * cb[1]; 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); 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) ) ); 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 information for all rings 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::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; //cout << ringShiftSum << " " << g << " " << Rings[i].resID << 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() { Mol conf = Conformers[1]; Mol::iterator it; //retreive the donor and acceptor list from coordinates for ( it = conf.begin(); it != conf.end(); it++ ) { if( isAcceptor(it->second) ) acceptorList.push_back(it->first); else { PDB_Entry temp = isDonor(it->second); if(temp.atomName.empty()) continue; donorList[it->first] = temp.atomNum; } } float dist, angle; map::iterator itD; for(int i = 0; i < acceptorList.size(); i++) //loop over acceptor list { PDB_Entry A = conf[ acceptorList[i] ]; for( itD = donorList.begin(); itD != donorList.end(); itD++ ) { PDB_Entry D = conf[ itD->first ], D_Heavy = conf[ itD->second ]; dist = getDist(D,A); angle = getBondAngle(D, D_Heavy, A); if( dist <= 3.5 && fabs(angle) <= 35) { HBDistList[A.resNum][A.atomName] = dist; HBDistList[D.resNum][D.atomName] = dist; } } } } // 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]; } bool PDB::isAcceptor(PDB_Entry A) { if(A.atomName == "O" || A.atomName == "N" || A.atomName == "OT1" || A.atomName == "OT2" ) return true; if(A.resName == "ASN" && A.atomName == "OD1") return true; if(A.resName == "ASP" && (A.atomName == "OD1" || A.atomName == "OD2") ) return true; if((A.resName == "CYS"||A.resName == "cys") && A.atomName == "SG") return true; if(A.resName == "GLN" && A.atomName == "OE1") return true; if(A.resName == "GLU" && (A.atomName == "OE1" || A.atomName == "OE2") ) return true; if(A.resName == "HIS" && (A.atomName == "ND1" || A.atomName == "NE1") ) return true; if(A.resName == "MET" && A.atomName == "SD") return true; if(A.resName == "SER" && A.atomName == "OG") return true; if(A.resName == "THR" && A.atomName == "OG1") return true; if(A.resName == "TYR" && A.atomName == "OH") return true; return false; } PDB_Entry PDB::isDonor(PDB_Entry D) { if(D.atomName == "HN") return ATOMS[1][D.resNum]["N"]; 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; } int PDB::getVectors(DoubleD2 &vec, int rr1, int rrN, vector &atomNames, vector &segList) { int k = 0; int n = segList.size()/2; for(int i = rr1; i <= rrN; i++) { bool flag = false; for(int j = 0; j < n; j++) { if( i>= segList[2*j] && i <= segList[2*j+1]) flag = true; } if(!flag && n>0) continue; //cout << "**** residue" << i << endl; for(int n = 0; n < atomNames.size(); n++) { if( residList[i].length() == 0 ) { cerr << "Missed residue " << i << "\t" << PDBfileName << endl; return k; } PDB_Entry a = getEntry(1,i,atomNames[n]); if( a.atomName == atomNames[n] ) k++; else { cerr << "Missed atom " << atomNames[n] << " @" << i << "\t" << PDBfileName << endl; return k; } vec[1][k] = a.X; vec[2][k] = a.Y; vec[3][k] = a.Z; } } return k; //cout << k << endl; } //get fragment from resID 1 to last residue int PDB::getVectors(DoubleD2 &vec, vector &atomNames, vector &segList) { return getVectors(vec, r1, rN, atomNames, segList); } double PDB::r_sign(double x, double y) { if( y >= 0 ) return fabs(x); else return -fabs(x); } // ================================================================== // SVDCMP: Singular value decomposition. // // Lit.: W. H. Press et. al., Numerical Recipes, // Cambridge University Press (1986), p. 60. // ================================================================== int PDB::svd_(DoubleD2 &a, DoubleD1 &w, DoubleD2 &v) { double cc, f, g, hh; int i, j, k, l; double s, x, y, zz; int jj, nm; double rv1[100]; int its; double scale, anorm; l = 0; g = 0.f; scale = 0.f; anorm = 0.f; for (i = 1; i <= 3; ++i) { l = i + 1; rv1[i - 1] = scale * g; g = 0.f; s = 0.f; scale = 0.f; if (i <= 3) { for (k = i; k <= 3; ++k) scale += fabs( a[k][i] ); if (scale != 0.f) { for (k = i; k <= 3; ++k) { a[k][i] /= scale; s += a[k][i] * a[k][i]; } f = a[i][i]; g = -r_sign(sqrt(s), f); hh = f * g - s; a[i][i] = f - g; if (i != 3) { for (j = l; j <= 3; ++j) { s = 0.f; for (k = i; k <= 3; ++k) s += a[k][i] * a[k][j]; f = s / hh; for (k = i; k <= 3; ++k) a[k][j] += f * a[k][i]; } } for (k = i; k <= 3; ++k) a[k][i] = scale * a[k][i]; } } w[i] = scale * g; g = 0.f; s = 0.f; scale = 0.f; if (i <= 3 && i != 3) { for (k = l; k <= 3; ++k) scale += fabs( a[i][k] ); if (scale != 0.f) { for (k = l; k <= 3; ++k) { a[i][k] /= scale; s += a[i][k] * a[i][k]; } f = a[i][l]; g = -r_sign(sqrt(s), f); hh = f * g - s; a[i][l] = f - g; for (k = l; k <= 3; ++k) rv1[k - 1] = a[i][k] / hh; if (i != 3) { for (j = l; j <= 3; ++j) { s = 0.f; for (k = l; k <= 3; ++k) s += a[j][k] * a[i][k]; for (k = l; k <= 3; ++k) a[j][k] += s * rv1[k - 1]; } } for (k = l; k <= 3; ++k) { a[i][k] = scale * a[i][k]; } } } anorm = max( anorm, fabs(w[i]) + fabs(rv1[i - 1]) ); } for (i = 3; i >= 1; --i) { if (i < 3) { if (g != 0.f) { for (j = l; j <= 3; ++j) v[j][i] = a[i][j] / a[i][l] / g; for (j = l; j <= 3; ++j) { s = 0.f; for (k = l; k <= 3; ++k) s += a[i][k] * v[k][j]; for (k = l; k <= 3; ++k) v[k][j] += s * v[k][i]; } } for (j = l; j <= 3; ++j) { v[i][j] = 0.f; v[j][i] = 0.f; } } v[i][i] = 1.f; g = rv1[i - 1]; l = i; } for (i = 3; i >= 1; --i) { l = i + 1; g = w[i]; if (i < 3) { for (j = l; j <= 3; ++j) a[i][j] = 0.f; } if (g != 0.f) { g = 1.f / g; if (i != 3) { for (j = l; j <= 3; ++j) { s = 0.f; for (k = l; k <= 3; ++k) s += a[k][i] * a[k][j]; f = s / a[i][i] * g; for (k = i; k <= 3; ++k) a[k][j] += f * a[k][i]; } } for (j = i; j <= 3; ++j) a[j][i] *= g; } else { for (j = i; j <= 3; ++j) a[j][i] = 0.f; } a[i][i] += 1.f; } for (k = 3; k >= 1; --k) { for (its = 1; its <= 30; ++its) { nm = 0; for (l = k; l >= 1; --l) { nm = l - 1; if ( (fabs(rv1[l - 1]) + anorm) == anorm) goto L2; if ( (fabs(w[nm]) + anorm) == anorm) goto L1; } L1: cc = 0.f; s = 1.f; for (i = l; i <= k; ++i) { f = s * rv1[i - 1]; if (fabs(f) + anorm != anorm) { g = w[i]; hh = sqrt(f * f + g * g); w[i] = hh; hh = 1.f / hh; cc = g * hh; s = -(f * hh); for (j = 1; j <= 3; ++j) { y = a[j][nm]; zz = a[j][i]; a[j][nm] = y * cc + zz * s; a[j][i] = -(y * s) + zz * cc; } } } L2: zz = w[k]; if (l == k) { if (zz < 0.f) { w[k] = -zz; for (j = 1; j <= 3; ++j) v[j][k] = -v[j][k]; } goto L3; } if (its == 30) printf("SVDCMP: No convergence in 30 iterations\n"); x = w[l]; nm = k - 1; y = w[nm]; g = rv1[nm - 1]; hh = rv1[k - 1]; f = ((y - zz) * (y + zz) + (g - hh) * (g + hh)) / (hh * 2.f * y); /* PG avoid overflow */ /* PG G=SQRT(F*F+1.0) */ /* Computing 2nd power */ g = f * sqrt((1.f / f) * (1.f / f) + 1.f); f = ((x - zz) * (x + zz) + hh * (y / (f + r_sign(g, f)) - hh)) / x; cc = 1.f; s = 1.f; for (j = l; j <= nm; ++j) { i = j + 1; g = rv1[i - 1]; y = w[i]; hh = s * g; g = cc * g; zz = sqrt(f * f + hh * hh); rv1[j - 1] = zz; cc = f / zz; s = hh / zz; f = x * cc + g * s; g = -(x * s) + g * cc; hh = y * s; y *= cc; for (jj = 1; jj <= 3; ++jj) { x = v[jj][j]; zz = v[jj][i]; v[jj][j] = x * cc + zz * s; v[jj][i] = -(x * s) + zz * cc; } zz = sqrt(f * f + hh * hh); w[j] = zz; if (zz != 0.f) { zz = 1.f / zz; cc = f * zz; s = hh * zz; } f = cc * g + s * y; x = -(s * g) + cc * y; for (jj = 1; jj <= 3; ++jj) { y = a[jj][j]; zz = a[jj][i]; a[jj][j] = y * cc + zz * s; a[jj][i] = -(y * s) + zz * cc; } } rv1[l - 1] = 0.f; rv1[k - 1] = f; w[k] = x; } L3: ; } return 0; } //calculate minimized rmsd betweeb coordinates a and b //a and b: (3,n) matrices double PDB::rmsd(DoubleD2 &a, DoubleD2 &b, int n) { DoubleD1 e; //eigen value of SVD DoubleD1 a_c, b_c; //center of sets of coordinates a and b DoubleD2 as, bs; //relative coordinates of a and b to the center DoubleD2 u, v ; DoubleD2 rot_m; //rotation matrix for coordinate a for minimized rmsd int i, j, k; double rmsd = 0.0, detu, sum, one = 1.0; for(i = 1; i <= 3; i++) if( a[i].size() != n || b[i].size() != n ) { cerr << " Atom numbers are not same " << endl; return -1; } //calculate the center and relative coordinates for a and b for(i = 1; i <= 3; i++) { a_c[i] = 0.0; b_c[i] = 0.0; for(j = 1; j <= n; j++) { a_c[i] = a_c[i] + a[i][j]; b_c[i] = b_c[i] + b[i][j]; } a_c[i] = a_c[i] / (double)n; b_c[i] = b_c[i] / (double)n; for(j = 1; j <= n; j++) { as[i][j] = a[i][j] - a_c[i]; bs[i][j] = b[i][j] - b_c[i]; } } //calculate the 'distance matrix" from coordinates a to b for(j = 1; j <= 3; j++) { for(i = 1; i <= 3; i++) u[i][j] = 0.0; for(k = 1; k <= n; k++) { u[1][j] = u[1][j] + as[1][k] * bs[j][k]; u[2][j] = u[2][j] + as[2][k] * bs[j][k]; u[3][j] = u[3][j] + as[3][k] * bs[j][k]; } for(i = 1; i <= 3; i++) u[i][j] = u[i][j] / (double)n; } detu = u[1][1] * u[2][2] * u[3][3] + u[1][2] * u[2][3] * u[3][1]+ u[1][3] * u[2][1] * u[3][2] - u[1][3] * u[2][2] * u[3][1]- u[1][1] * u[2][3] * u[3][2] - u[1][2] * u[2][1] * u[3][3]; /* printf("ORIGIN:\n"); for(i = 1; i <= 3; i++){ for(j = 1; j <= 3; j++) printf("%6.2f ", u[i][j] ); printf("\n"); } */ //calculate the SVD for u svd_(u,e,v); /* printf("u:\n"); for(i = 1 ; i <= 3; i++){ for(j = 1; j <= 3; j++) printf("%6.2f ", u[i][j] ); printf("\n"); } printf("e:\n"); for(i = 1; i <= 3; i++) printf("%6.2f ", e[i] ); printf("\n"); printf("v:\n"); for(i = 1; i <= 3; i++) { for(j = 1; j <= 3; j++) printf("%6.2f ", v[i][j] ); printf("\n"); } */ sum = min(min(e[1],e[2]), e[3]); for(i = 1; i <= 3; i++) { if(e[i] == sum) e[i] = detu; e[i] = r_sign(1.0, e[i]); } //calculate the rotation matrix from coordinate a to b for minimized rmsd for(i = 1; i <= 3; i++) for(j = 1; j <= 3; j++) rot_m[i][j] = e[1]*v[i][1]*u[j][1] + e[2]*v[i][2]*u[j][2] + e[3]*v[i][3]*u[j][3]; //apply the rotation and calculate rmsd for(k = 1; k <= n; k++) { sum = rot_m[1][1] * as[1][k] + rot_m[1][2] * as[2][k] + rot_m[1][3] * as[3][k] - bs[1][k]; rmsd = rmsd + sum * sum; sum = rot_m[2][1] * as[1][k] + rot_m[2][2] * as[2][k] + rot_m[2][3] * as[3][k] - bs[2][k]; rmsd = rmsd + sum * sum; sum = rot_m[3][1] * as[1][k] + rot_m[3][2] * as[2][k] + rot_m[3][3 ]* as[3][k] - bs[3][k]; rmsd = rmsd + sum * sum; } return sqrt(rmsd/(double)n); } //calculate minimized rmsd betweeb coordinates a and b //a and b: (n,3) matrices double PDB::rmsd(vector &a, vector &b) { DoubleD2 as, bs; if( a.size() != b.size() ) return -1; if( a.size() == 0 ) return 0.0; // convert (n,3) matrices, which are indexed starting from 0, // to (3,n) matrices, which indexed starting from 1 for(int n = 0; n < a.size(); n++) for(int i = 0; i < 3; i++) { as[i+1][n+1] = a[n][i]; bs[i+1][n+1] = b[n][i]; } return rmsd(as, bs, a.size()); }