/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /* * PdbParser.cpp * * Created on: Aug 17, 2010 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "Helpers/MemDebug.hpp" #include "Helpers/Assert.hpp" #include "Helpers/Log.hpp" #include "Helpers/Verbose.hpp" #include "PdbParser.hpp" #include "World.hpp" #include "atom.hpp" #include "bond.hpp" #include "element.hpp" #include "molecule.hpp" #include "periodentafel.hpp" #include "Descriptors/AtomIdDescriptor.hpp" #include #include #include #include using namespace std; /** * Constructor. */ PdbParser::PdbParser() { knownKeys[" "] = PdbKey::noKey; // with this we can detect invalid keys knownKeys["x"] = PdbKey::x; knownKeys["Id"] = PdbKey::Id; knownKeys["Type"] = PdbKey::Type; knownKeys["extType"] = PdbKey::extType; knownKeys["name"] = PdbKey::name; knownKeys["resName"] = PdbKey::resName; knownKeys["chainID"] = PdbKey::chainID; knownKeys["resSeq"] = PdbKey::resSeq; knownKeys["occupancy"] = PdbKey::occupancy; knownKeys["tempFactor"] = PdbKey::tempFactor; knownKeys["segID"] = PdbKey::segID; knownKeys["charge"] = PdbKey::charge; } /** * Destructor. */ PdbParser::~PdbParser() { additionalAtomData.clear(); atomIdMap.clear(); knownKeys.clear(); } /** * Loads atoms from a tremolo-formatted file. * * \param tremolo file */ void PdbParser::load(istream* file) { // TODO: PdbParser::load implementation ASSERT(false, "Not implemented yet"); // string line; // string::size_type location; // // usedFields.clear(); // while (file->good()) { // std::getline(*file, line, '\n'); // if (usedFields.empty()) { // location = line.find("ATOMDATA", 0); // if (location != string::npos) { // parseAtomDataKeysLine(line, location + 8); // } // } // if (line.length() > 0 && line.at(0) != '#') { // readAtomDataLine(line); // } // } // // processNeighborInformation(); // adaptImprData(); // adaptTorsion(); } /** * Saves the World's current state into as a tremolo file. * * \param file where to save the state */ void PdbParser::save(ostream* file) { DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl); { // add initial remark *file << "REMARK created by molecuilder on "; time_t now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time' // ctime ends in \n\0, we have to cut away the newline std::string time(ctime(&now)); size_t pos = time.find('\n'); if (pos != 0) *file << time.substr(0,pos); else *file << time; *file << endl; } { vector AtomList = World::getInstance().getAllAtoms(); std::vector elementNo(MAX_ELEMENTS,1); char name[MAXSTRINGSIZE]; // write ATOMs int AtomNo = 1; // serial number starts at 1 in pdb int MolNo = 1; // residue number starts at 1 in pdb for (vector::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) { const size_t Z = (*atomIt)->getType()->getAtomicNumber(); sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), elementNo[Z]); elementNo[Z] = (elementNo[Z]+1) % 100; // confine to two digits const molecule *mol = (*atomIt)->getMolecule(); if (mol == NULL) { // for homeless atoms, MolNo = -1 is reserved MolNo = -1; } else { MolNo = mol->getId(); } saveLine(file, *atomIt, name, AtomNo, MolNo); setAtomId((*atomIt)->getId(), AtomNo); AtomNo++; } // write CONECTs for (vector::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) { writeNeighbors(file, 4, *atomIt); } } // END *file << "END" << endl; } /** * Writes one line of tremolo-formatted data to the provided stream. * * \param stream where to write the line to * \param *currentAtom the atom of which information should be written * \param *name name of atom, i.e. H01 * \param AtomNo serial number of atom * \param ResidueNo number of residue */ void PdbParser::saveLine(ostream* file, const atom* currentAtom, const char *name, const int AtomNo, const int ResidueNo) { *file << "ATOM "; *file << setw(6) << AtomNo; /* atom serial number */ *file << setw(1) << " "; *file << setfill(' ') << left << setw(4) << name << right; /* atom name */ *file << setw(1) << " "; *file << setfill(' ') << setw(3) << ((currentAtom->getMolecule() != NULL) ? currentAtom->getMolecule()->getName().substr(0,3) : "-"); /* residue name */ *file << setw(1) << " "; *file << setfill(' ') << setw(1) << (char)('a'+(unsigned char)(AtomNo % 26)); /* letter for chain */ *file << setw(4) << ResidueNo; /* residue sequence number */ *file << setw(4) << " "; for (int i=0;iat(i); /* positional coordinate in Angstroem */ } *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getValence(); /* occupancy */ *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getNoValenceOrbitals(); /* temperature factor */ *file << noshowpoint; *file << setw(6) << " "; *file << setw(4) << "0"; *file << setfill(' ') << setw(2) << currentAtom->getType()->getSymbol(); *file << setw(2) << "0"; *file << endl; } /** * Writes the neighbor information of one atom to the provided stream. * * \param *file where to write neighbor information to * \param MaxnumberOfNeighbors of neighbors * \param *currentAtom to the atom of which to take the neighbor information */ void PdbParser::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) { if (!currentAtom->ListOfBonds.empty()) { *file << "CONECT"; *file << setw(5) << getAtomId(currentAtom->getId()); int MaxNo = 0; for(BondList::iterator currentBond = currentAtom->ListOfBonds.begin(); currentBond != currentAtom->ListOfBonds.end(); ++currentBond) { if (MaxNo < MaxnumberOfNeighbors) { *file << setw(5) << getAtomId((*currentBond)->GetOtherAtom(currentAtom)->getId()); } MaxNo++; } *file << endl; } } /** Retrieves a value from PdbParser::atomIdMap. * \param atomid key * \return value */ int PdbParser::getAtomId(int atomid) const { ASSERT(atomIdMap.find(atomid) != atomIdMap.end(), "PdbParser::getAtomId: atomid not present in Map."); return (atomIdMap.find(atomid)->second); } /** Sets an entry in PdbParser::atomIdMap. * \param localatomid key * \param atomid value * \return true - key not present, false - value present */ void PdbParser::setAtomId(int localatomid, int atomid) { pair::iterator, bool > inserter; inserter = atomIdMap.insert( pair(localatomid, atomid) ); ASSERT(inserter.second, "PdbParser::setAtomId: atomId already present in Map."); } /** * Reads one data line of a tremolo file and interprets it according to the keys * obtained from the ATOMDATA line. * * \param line to parse as an atom */ void PdbParser::readAtomDataLine(string line) { // vector::iterator it; // stringstream lineStream; // atom* newAtom = World::getInstance().createAtom(); // PdbAtomInfoContainer *atomInfo = NULL; // additionalAtomData[newAtom->getId()] = *(new PdbAtomInfoContainer); // atomInfo = &additionalAtomData[newAtom->getId()]; // PdbKey::atomDataKey currentField; // string word; // int oldId; // double tmp; // // lineStream << line; // for (it = usedFields.begin(); it < usedFields.end(); it++) { // currentField = knownKeys[it->substr(0, it->find("="))]; // switch (currentField) { // case PdbKey::x : // // for the moment, assume there are always three dimensions // for (int i=0;i> tmp; // newAtom->set(i, tmp); // } // break; // case PdbKey::u : // // for the moment, assume there are always three dimensions // lineStream >> newAtom->AtomicVelocity[0]; // lineStream >> newAtom->AtomicVelocity[1]; // lineStream >> newAtom->AtomicVelocity[2]; // break; // case PdbKey::Type : // char type[3]; // lineStream >> type; // newAtom->setType(World::getInstance().getPeriode()->FindElement(type)); // ASSERT(newAtom->getType(), "Type was not set for this atom"); // break; // case PdbKey::Id : // lineStream >> oldId; // atomIdMap[oldId] = newAtom->getId(); // break; // case PdbKey::neighbors : // readNeighbors(&lineStream, // atoi(it->substr(it->find("=") + 1, 1).c_str()), newAtom->getId()); // break; // default : // lineStream >> word; // atomInfo->set(currentField, word); // break; // } // } } /** * Reads neighbor information for one atom from the input. * * \param stream where to read the information from * \param number of neighbors to read * \param world id of the atom the information belongs to */ void PdbParser::readNeighbors(stringstream* line, int numberOfNeighbors, int atomId) { // int neighborId = 0; // for (int i = 0; i < numberOfNeighbors; i++) { // *line >> neighborId; // // 0 is used to fill empty neighbor positions in the tremolo file. // if (neighborId > 0) { // additionalAtomData[atomId].neighbors.push_back(neighborId); // } // } } /** * Adds the collected neighbor information to the atoms in the world. The atoms * are found by their current ID and mapped to the corresponding atoms with the * Id found in the parsed file. */ void PdbParser::processNeighborInformation() { // if (!isUsedField("neighbors")) { // return; // } // // for(map::iterator currentInfo = additionalAtomData.begin(); // currentInfo != additionalAtomData.end(); currentInfo++ // ) { // for(vector::iterator neighbor = currentInfo->second.neighbors.begin(); // neighbor != currentInfo->second.neighbors.end(); neighbor++ // ) { // World::getInstance().getAtom(AtomById(currentInfo->first)) // ->addBond(World::getInstance().getAtom(AtomById(atomIdMap[*neighbor]))); // } // } } /** * Replaces atom IDs read from the file by the corresponding world IDs. All IDs * IDs of the input string will be replaced; expected separating characters are * "-" and ",". * * \param string in which atom IDs should be adapted * * \return input string with modified atom IDs */ string PdbParser::adaptIdDependentDataString(string data) { // // there might be no IDs // if (data == "-") { // return "-"; // } // // char separator; // int id; // stringstream line, result; // line << data; // // line >> id; // result << atomIdMap[id]; // while (line.good()) { // line >> separator >> id; // result << separator << atomIdMap[id]; // } // // return result.str(); return ""; } PdbAtomInfoContainer::PdbAtomInfoContainer() : name("-"), resName("-"), chainID("0"), resSeq("0"), occupancy("0"), tempFactor("0"), segID("0"), charge("0") {} void PdbAtomInfoContainer::set(PdbKey::PdbDataKey key, string value) { switch (key) { case PdbKey::extType : extType = value; break; case PdbKey::name : name = value; break; case PdbKey::resName : resName = value; break; case PdbKey::chainID : chainID = value; break; case PdbKey::resSeq : resSeq = value; break; case PdbKey::occupancy : occupancy = value; break; case PdbKey::tempFactor : tempFactor = value; break; case PdbKey::segID : segID = value; break; case PdbKey::charge : charge = value; break; default : cout << "Unknown key: " << key << ", value: " << value << endl; break; } } string PdbAtomInfoContainer::get(PdbKey::PdbDataKey key) { switch (key) { case PdbKey::extType : return extType; case PdbKey::name : return name; case PdbKey::resName : return resName; case PdbKey::chainID : return chainID; case PdbKey::resSeq : return resSeq; case PdbKey::occupancy : return occupancy; case PdbKey::tempFactor : return tempFactor; case PdbKey::segID : return segID; case PdbKey::charge : return charge; default : cout << "Unknown key: " << key << endl; return ""; } }