[3ae731] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
---|
| 5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | */
|
---|
| 7 |
|
---|
| 8 | /*
|
---|
| 9 | * PdbParser.cpp
|
---|
| 10 | *
|
---|
| 11 | * Created on: Aug 17, 2010
|
---|
| 12 | * Author: heber
|
---|
| 13 | */
|
---|
| 14 |
|
---|
| 15 | // include config.h
|
---|
| 16 | #ifdef HAVE_CONFIG_H
|
---|
| 17 | #include <config.h>
|
---|
| 18 | #endif
|
---|
| 19 |
|
---|
[ad011c] | 20 | #include "CodePatterns/MemDebug.hpp"
|
---|
[3ae731] | 21 |
|
---|
[ad011c] | 22 | #include "CodePatterns/Assert.hpp"
|
---|
| 23 | #include "CodePatterns/Log.hpp"
|
---|
| 24 | #include "CodePatterns/toString.hpp"
|
---|
| 25 | #include "CodePatterns/Verbose.hpp"
|
---|
[3ae731] | 26 | #include "atom.hpp"
|
---|
| 27 | #include "bond.hpp"
|
---|
[bb6193] | 28 | #include "element.hpp"
|
---|
| 29 | #include "molecule.hpp"
|
---|
[3ae731] | 30 | #include "periodentafel.hpp"
|
---|
| 31 | #include "Descriptors/AtomIdDescriptor.hpp"
|
---|
[4fbca9c] | 32 | #include "Parser/PdbParser.hpp"
|
---|
[073a9e4] | 33 | #include "World.hpp"
|
---|
| 34 | #include "WorldTime.hpp"
|
---|
[bb6193] | 35 |
|
---|
[3ae731] | 36 | #include <map>
|
---|
| 37 | #include <vector>
|
---|
| 38 |
|
---|
[bb6193] | 39 | #include <iostream>
|
---|
| 40 | #include <iomanip>
|
---|
[3ae731] | 41 |
|
---|
| 42 | using namespace std;
|
---|
| 43 |
|
---|
| 44 | /**
|
---|
| 45 | * Constructor.
|
---|
| 46 | */
|
---|
| 47 | PdbParser::PdbParser() {
|
---|
[4fbca9c] | 48 | knownTokens["ATOM"] = PdbKey::Atom;
|
---|
[16462f] | 49 | knownTokens["HETATM"] = PdbKey::Atom;
|
---|
[4fbca9c] | 50 | knownTokens["TER"] = PdbKey::Filler;
|
---|
| 51 | knownTokens["END"] = PdbKey::EndOfFile;
|
---|
| 52 | knownTokens["CONECT"] = PdbKey::Connect;
|
---|
| 53 | knownTokens["REMARK"] = PdbKey::Remark;
|
---|
| 54 | knownTokens[""] = PdbKey::EndOfFile;
|
---|
[16462f] | 55 |
|
---|
| 56 | // argh, why can't just PdbKey::X+(size_t)i
|
---|
| 57 | PositionEnumMap[0] = PdbKey::X;
|
---|
| 58 | PositionEnumMap[1] = PdbKey::Y;
|
---|
| 59 | PositionEnumMap[2] = PdbKey::Z;
|
---|
[3ae731] | 60 | }
|
---|
| 61 |
|
---|
| 62 | /**
|
---|
| 63 | * Destructor.
|
---|
| 64 | */
|
---|
| 65 | PdbParser::~PdbParser() {
|
---|
| 66 | additionalAtomData.clear();
|
---|
| 67 | atomIdMap.clear();
|
---|
[4fbca9c] | 68 | }
|
---|
| 69 |
|
---|
| 70 |
|
---|
| 71 | /** Parses the initial word of the given \a line and returns the token type.
|
---|
| 72 | *
|
---|
| 73 | * @param line line to scan
|
---|
| 74 | * @return token type
|
---|
| 75 | */
|
---|
| 76 | enum PdbKey::KnownTokens PdbParser::getToken(string &line)
|
---|
| 77 | {
|
---|
| 78 | // look for first space
|
---|
| 79 | const size_t space_location = line.find(' ');
|
---|
| 80 | const size_t tab_location = line.find('\t');
|
---|
| 81 | size_t location = space_location < tab_location ? space_location : tab_location;
|
---|
| 82 | string token;
|
---|
| 83 | if (location != string::npos) {
|
---|
| 84 | //DoLog(1) && (Log() << Verbose(1) << "Found space at position " << space_location << std::endl);
|
---|
| 85 | token = line.substr(0,space_location);
|
---|
| 86 | } else {
|
---|
| 87 | token = line;
|
---|
| 88 | }
|
---|
| 89 |
|
---|
| 90 | //DoLog(1) && (Log() << Verbose(1) << "Token is " << token << std::endl);
|
---|
| 91 | if (knownTokens.count(token) == 0)
|
---|
| 92 | return PdbKey::NoToken;
|
---|
| 93 | else
|
---|
| 94 | return knownTokens[token];
|
---|
| 95 |
|
---|
| 96 | return PdbKey::NoToken;
|
---|
[3ae731] | 97 | }
|
---|
| 98 |
|
---|
| 99 | /**
|
---|
[4fbca9c] | 100 | * Loads atoms from a PDB-formatted file.
|
---|
[3ae731] | 101 | *
|
---|
[4fbca9c] | 102 | * \param PDB file
|
---|
[3ae731] | 103 | */
|
---|
| 104 | void PdbParser::load(istream* file) {
|
---|
[4fbca9c] | 105 | string line;
|
---|
| 106 | size_t linecount = 0;
|
---|
| 107 | enum PdbKey::KnownTokens token;
|
---|
| 108 |
|
---|
[16462f] | 109 | // reset atomIdMap for this file (to correctly parse CONECT entries)
|
---|
| 110 | atomIdMap.clear();
|
---|
| 111 |
|
---|
[4fbca9c] | 112 | molecule *newmol = World::getInstance().createMolecule();
|
---|
| 113 | newmol->ActiveFlag = true;
|
---|
| 114 | bool NotEndOfFile = true;
|
---|
| 115 | // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
|
---|
| 116 | World::getInstance().getMolecules()->insert(newmol);
|
---|
| 117 | while (NotEndOfFile) {
|
---|
| 118 | std::getline(*file, line, '\n');
|
---|
| 119 | // extract first token
|
---|
| 120 | token = getToken(line);
|
---|
| 121 | //DoLog(1) && (Log() << Verbose(1) << " Recognized token of type : " << token << std::endl);
|
---|
| 122 | switch (token) {
|
---|
| 123 | case PdbKey::Atom:
|
---|
| 124 | readAtomDataLine(line, newmol);
|
---|
| 125 | break;
|
---|
| 126 | case PdbKey::Remark:
|
---|
| 127 | break;
|
---|
| 128 | case PdbKey::Connect:
|
---|
| 129 | readNeighbors(line);
|
---|
| 130 | break;
|
---|
| 131 | case PdbKey::Filler:
|
---|
| 132 | break;
|
---|
| 133 | case PdbKey::EndOfFile:
|
---|
| 134 | NotEndOfFile = false;
|
---|
| 135 | break;
|
---|
| 136 | default:
|
---|
| 137 | // TODO: put a throw here
|
---|
| 138 | DoeLog(2) && (eLog() << Verbose(2) << "Unknown token: '" << line << "'" << std::endl);
|
---|
| 139 | //ASSERT(0, "PdbParser::load() - Unknown token in line "+toString(linecount)+": "+line+".");
|
---|
| 140 | break;
|
---|
| 141 | }
|
---|
| 142 | NotEndOfFile = NotEndOfFile && (file->good());
|
---|
| 143 | linecount++;
|
---|
| 144 | }
|
---|
[3ae731] | 145 | }
|
---|
| 146 |
|
---|
| 147 | /**
|
---|
[73916f] | 148 | * Saves the \a atoms into as a PDB file.
|
---|
[3ae731] | 149 | *
|
---|
| 150 | * \param file where to save the state
|
---|
[73916f] | 151 | * \param atoms atoms to store
|
---|
[3ae731] | 152 | */
|
---|
[73916f] | 153 | void PdbParser::save(ostream* file, const std::vector<atom *> &AtomList)
|
---|
| 154 | {
|
---|
[bb6193] | 155 | DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl);
|
---|
| 156 | {
|
---|
| 157 | // add initial remark
|
---|
| 158 | *file << "REMARK created by molecuilder on ";
|
---|
| 159 | time_t now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time'
|
---|
| 160 | // ctime ends in \n\0, we have to cut away the newline
|
---|
| 161 | std::string time(ctime(&now));
|
---|
| 162 | size_t pos = time.find('\n');
|
---|
| 163 | if (pos != 0)
|
---|
| 164 | *file << time.substr(0,pos);
|
---|
| 165 | else
|
---|
| 166 | *file << time;
|
---|
| 167 | *file << endl;
|
---|
| 168 | }
|
---|
[3ae731] | 169 |
|
---|
[16462f] | 170 | // we distribute serials, hence clear map beforehand
|
---|
[4fbca9c] | 171 | atomIdMap.clear();
|
---|
[bb6193] | 172 | {
|
---|
[16462f] | 173 | std::map<size_t,size_t> MolIdMap;
|
---|
| 174 | size_t MolNo = 1; // residue number starts at 1 in pdb
|
---|
[73916f] | 175 | for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
|
---|
| 176 | const molecule *mol = (*atomIt)->getMolecule();
|
---|
| 177 | if ((mol != NULL) && (MolIdMap.find(mol->getId()) == MolIdMap.end())) {
|
---|
| 178 | MolIdMap[mol->getId()] = MolNo++;
|
---|
| 179 | }
|
---|
[16462f] | 180 | }
|
---|
| 181 | const size_t MaxMol = MolNo;
|
---|
[bb6193] | 182 |
|
---|
[16462f] | 183 | // have a count per element and per molecule (0 is for all homeless atoms)
|
---|
| 184 | std::vector<int> **elementNo = new std::vector<int>*[MaxMol];
|
---|
| 185 | for (size_t i = 0; i < MaxMol; ++i)
|
---|
| 186 | elementNo[i] = new std::vector<int>(MAX_ELEMENTS,1);
|
---|
[bb6193] | 187 | char name[MAXSTRINGSIZE];
|
---|
[16462f] | 188 | std::string ResidueName;
|
---|
[bb6193] | 189 |
|
---|
| 190 | // write ATOMs
|
---|
| 191 | int AtomNo = 1; // serial number starts at 1 in pdb
|
---|
[73916f] | 192 | for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
|
---|
[93fd43e] | 193 | PdbAtomInfoContainer &atomInfo = getadditionalAtomData(*atomIt);
|
---|
[16462f] | 194 | // gather info about residue
|
---|
[bb6193] | 195 | const molecule *mol = (*atomIt)->getMolecule();
|
---|
[16462f] | 196 | if (mol == NULL) {
|
---|
| 197 | MolNo = 0;
|
---|
| 198 | atomInfo.set(PdbKey::resSeq, "0");
|
---|
[bb6193] | 199 | } else {
|
---|
[16462f] | 200 | ASSERT(MolIdMap.find(mol->getId()) != MolIdMap.end(),
|
---|
| 201 | "PdbParser::save() - Mol id "+toString(mol->getId())+" not present despite we set it?!");
|
---|
| 202 | MolNo = MolIdMap[mol->getId()];
|
---|
| 203 | atomInfo.set(PdbKey::resSeq, toString(MolIdMap[mol->getId()]));
|
---|
| 204 | if (atomInfo.get<std::string>(PdbKey::resName) == "-")
|
---|
| 205 | atomInfo.set(PdbKey::resName, mol->getName().substr(0,3));
|
---|
| 206 | }
|
---|
| 207 | // get info about atom
|
---|
| 208 | const size_t Z = (*atomIt)->getType()->getAtomicNumber();
|
---|
| 209 | if (atomInfo.get<std::string>(PdbKey::name) == "-") { // if no name set, give it a new name
|
---|
| 210 | sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), (*elementNo[MolNo])[Z]);
|
---|
| 211 | (*elementNo[MolNo])[Z] = ((*elementNo[MolNo])[Z]+1) % 100; // confine to two digits
|
---|
| 212 | atomInfo.set(PdbKey::name, name);
|
---|
[bb6193] | 213 | }
|
---|
[16462f] | 214 | // set position
|
---|
| 215 | for (size_t i=0; i<NDIM;++i) {
|
---|
| 216 | stringstream position;
|
---|
| 217 | position << setw(8) << fixed << setprecision(3) << (*atomIt)->getPosition().at(i);
|
---|
| 218 | atomInfo.set(PositionEnumMap[i], position.str());
|
---|
| 219 | }
|
---|
| 220 | // change element and charge if changed
|
---|
| 221 | if (atomInfo.get<std::string>(PdbKey::element) != (*atomIt)->getType()->getSymbol())
|
---|
| 222 | atomInfo.set(PdbKey::element, (*atomIt)->getType()->getSymbol());
|
---|
| 223 | setSerial((*atomIt)->getId(), AtomNo);
|
---|
| 224 | atomInfo.set(PdbKey::serial, toString(AtomNo));
|
---|
| 225 |
|
---|
| 226 | // finally save the line
|
---|
| 227 | saveLine(file, atomInfo);
|
---|
[bb6193] | 228 | AtomNo++;
|
---|
| 229 | }
|
---|
[16462f] | 230 | for (size_t i = 0; i < MaxMol; ++i)
|
---|
| 231 | delete elementNo[i];
|
---|
| 232 | delete elementNo;
|
---|
[3ae731] | 233 |
|
---|
[bb6193] | 234 | // write CONECTs
|
---|
[73916f] | 235 | for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
|
---|
[bb6193] | 236 | writeNeighbors(file, 4, *atomIt);
|
---|
| 237 | }
|
---|
[3ae731] | 238 | }
|
---|
| 239 |
|
---|
[bb6193] | 240 | // END
|
---|
| 241 | *file << "END" << endl;
|
---|
[3ae731] | 242 | }
|
---|
| 243 |
|
---|
[93fd43e] | 244 | /** Either returns reference to present entry or creates new with default values.
|
---|
| 245 | *
|
---|
| 246 | * @param _atom atom whose entry we desire
|
---|
| 247 | * @return
|
---|
| 248 | */
|
---|
| 249 | PdbAtomInfoContainer& PdbParser::getadditionalAtomData(atom *_atom)
|
---|
| 250 | {
|
---|
| 251 | if (additionalAtomData.find(_atom->getId()) != additionalAtomData.end()) {
|
---|
| 252 | } else if (additionalAtomData.find(_atom->father->getId()) != additionalAtomData.end()) {
|
---|
| 253 | // use info from direct father
|
---|
| 254 | additionalAtomData[_atom->getId()] = additionalAtomData[_atom->father->getId()];
|
---|
| 255 | } else if (additionalAtomData.find(_atom->GetTrueFather()->getId()) != additionalAtomData.end()) {
|
---|
| 256 | // use info from topmost father
|
---|
| 257 | additionalAtomData[_atom->getId()] = additionalAtomData[_atom->GetTrueFather()->getId()];
|
---|
| 258 | } else {
|
---|
| 259 | // create new entry use default values if nothing else is known
|
---|
| 260 | additionalAtomData[_atom->getId()] = defaultAdditionalData;
|
---|
| 261 | }
|
---|
| 262 | return additionalAtomData[_atom->getId()];
|
---|
| 263 | }
|
---|
| 264 |
|
---|
[3ae731] | 265 | /**
|
---|
[4fbca9c] | 266 | * Writes one line of PDB-formatted data to the provided stream.
|
---|
[3ae731] | 267 | *
|
---|
| 268 | * \param stream where to write the line to
|
---|
[bb6193] | 269 | * \param *currentAtom the atom of which information should be written
|
---|
| 270 | * \param AtomNo serial number of atom
|
---|
[16462f] | 271 | * \param *name name of atom, i.e. H01
|
---|
| 272 | * \param ResidueName Name of molecule
|
---|
[bb6193] | 273 | * \param ResidueNo number of residue
|
---|
[3ae731] | 274 | */
|
---|
[16462f] | 275 | void PdbParser::saveLine(
|
---|
| 276 | ostream* file,
|
---|
| 277 | const PdbAtomInfoContainer &atomInfo)
|
---|
| 278 | {
|
---|
| 279 | *file << setfill(' ') << left << setw(6)
|
---|
| 280 | << atomInfo.get<std::string>(PdbKey::token);
|
---|
| 281 | *file << setfill(' ') << right << setw(5)
|
---|
| 282 | << atomInfo.get<int>(PdbKey::serial); /* atom serial number */
|
---|
| 283 | *file << " "; /* char 12 is empty */
|
---|
| 284 | *file << setfill(' ') << left << setw(4)
|
---|
| 285 | << atomInfo.get<std::string>(PdbKey::name); /* atom name */
|
---|
| 286 | *file << setfill(' ') << left << setw(1)
|
---|
| 287 | << atomInfo.get<std::string>(PdbKey::altLoc); /* alternate location/conformation */
|
---|
| 288 | *file << setfill(' ') << left << setw(3)
|
---|
| 289 | << atomInfo.get<std::string>(PdbKey::resName); /* residue name */
|
---|
| 290 | *file << " "; /* char 21 is empty */
|
---|
| 291 | *file << setfill(' ') << left << setw(1)
|
---|
| 292 | << atomInfo.get<std::string>(PdbKey::chainID); /* chain identifier */
|
---|
| 293 | *file << setfill(' ') << left << setw(4)
|
---|
| 294 | << atomInfo.get<int>(PdbKey::resSeq); /* residue sequence number */
|
---|
| 295 | *file << setfill(' ') << left << setw(1)
|
---|
| 296 | << atomInfo.get<std::string>(PdbKey::iCode); /* iCode */
|
---|
| 297 | *file << " "; /* char 28-30 are empty */
|
---|
| 298 | // have the following operate on stringstreams such that format specifiers
|
---|
| 299 | // only act on these
|
---|
| 300 | for (size_t i=0;i<NDIM;++i) {
|
---|
| 301 | stringstream position;
|
---|
| 302 | position << fixed << setprecision(3) << showpoint
|
---|
| 303 | << atomInfo.get<double>(PositionEnumMap[i]);
|
---|
| 304 | *file << setfill(' ') << right << setw(8) << position.str();
|
---|
| 305 | }
|
---|
| 306 | {
|
---|
| 307 | stringstream occupancy;
|
---|
| 308 | occupancy << fixed << setprecision(2) << showpoint
|
---|
| 309 | << atomInfo.get<double>(PdbKey::occupancy); /* occupancy */
|
---|
| 310 | *file << setfill(' ') << right << setw(6) << occupancy.str();
|
---|
[3ae731] | 311 | }
|
---|
[16462f] | 312 | {
|
---|
| 313 | stringstream tempFactor;
|
---|
| 314 | tempFactor << fixed << setprecision(2) << showpoint
|
---|
| 315 | << atomInfo.get<double>(PdbKey::tempFactor); /* temperature factor */
|
---|
| 316 | *file << setfill(' ') << right << setw(6) << tempFactor.str();
|
---|
| 317 | }
|
---|
| 318 | *file << " "; /* char 68-76 are empty */
|
---|
| 319 | *file << setfill(' ') << right << setw(2) << atomInfo.get<std::string>(PdbKey::element); /* element */
|
---|
| 320 | *file << setfill(' ') << right << setw(2) << atomInfo.get<int>(PdbKey::charge); /* charge */
|
---|
[3ae731] | 321 |
|
---|
| 322 | *file << endl;
|
---|
| 323 | }
|
---|
| 324 |
|
---|
| 325 | /**
|
---|
| 326 | * Writes the neighbor information of one atom to the provided stream.
|
---|
| 327 | *
|
---|
[9d83b6] | 328 | * Note that ListOfBonds of WorldTime::CurrentTime is used.
|
---|
| 329 | *
|
---|
[bb6193] | 330 | * \param *file where to write neighbor information to
|
---|
| 331 | * \param MaxnumberOfNeighbors of neighbors
|
---|
| 332 | * \param *currentAtom to the atom of which to take the neighbor information
|
---|
[3ae731] | 333 | */
|
---|
[bb6193] | 334 | void PdbParser::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) {
|
---|
[4c1230] | 335 | int MaxNo = MaxnumberOfNeighbors;
|
---|
[9d83b6] | 336 | const BondList & ListOfBonds = currentAtom->getListOfBonds();
|
---|
| 337 | if (!ListOfBonds.empty()) {
|
---|
| 338 | for(BondList::const_iterator currentBond = ListOfBonds.begin(); currentBond != ListOfBonds.end(); ++currentBond) {
|
---|
[4c1230] | 339 | if (MaxNo >= MaxnumberOfNeighbors) {
|
---|
| 340 | *file << "CONECT";
|
---|
[16462f] | 341 | *file << setw(5) << getSerial(currentAtom->getId());
|
---|
[4c1230] | 342 | MaxNo = 0;
|
---|
[bb6193] | 343 | }
|
---|
[16462f] | 344 | *file << setw(5) << getSerial((*currentBond)->GetOtherAtom(currentAtom)->getId());
|
---|
[bb6193] | 345 | MaxNo++;
|
---|
[4c1230] | 346 | if (MaxNo == MaxnumberOfNeighbors)
|
---|
| 347 | *file << "\n";
|
---|
[3ae731] | 348 | }
|
---|
[4c1230] | 349 | if (MaxNo != MaxnumberOfNeighbors)
|
---|
| 350 | *file << "\n";
|
---|
[3ae731] | 351 | }
|
---|
| 352 | }
|
---|
| 353 |
|
---|
[21585f] | 354 | /** Retrieves a value from PdbParser::atomIdMap.
|
---|
| 355 | * \param atomid key
|
---|
| 356 | * \return value
|
---|
| 357 | */
|
---|
[4fbca9c] | 358 | size_t PdbParser::getSerial(const size_t atomid) const
|
---|
[21585f] | 359 | {
|
---|
| 360 | ASSERT(atomIdMap.find(atomid) != atomIdMap.end(), "PdbParser::getAtomId: atomid not present in Map.");
|
---|
| 361 | return (atomIdMap.find(atomid)->second);
|
---|
| 362 | }
|
---|
| 363 |
|
---|
| 364 | /** Sets an entry in PdbParser::atomIdMap.
|
---|
| 365 | * \param localatomid key
|
---|
| 366 | * \param atomid value
|
---|
| 367 | * \return true - key not present, false - value present
|
---|
| 368 | */
|
---|
[16462f] | 369 | void PdbParser::setSerial(const size_t localatomid, const size_t atomid)
|
---|
[21585f] | 370 | {
|
---|
[4fbca9c] | 371 | pair<std::map<size_t,size_t>::iterator, bool > inserter;
|
---|
[16462f] | 372 | // DoLog(1) && (Log() << Verbose(1) << "PdbParser::setAtomId() - Inserting ("
|
---|
| 373 | // << localatomid << " -> " << atomid << ")." << std::endl);
|
---|
[4fbca9c] | 374 | inserter = atomIdMap.insert( make_pair(localatomid, atomid) );
|
---|
[21585f] | 375 | ASSERT(inserter.second, "PdbParser::setAtomId: atomId already present in Map.");
|
---|
| 376 | }
|
---|
| 377 |
|
---|
[4fbca9c] | 378 | /** Parse an ATOM line from a PDB file.
|
---|
| 379 | *
|
---|
| 380 | * Reads one data line of a pdstatus file and interprets it according to the
|
---|
| 381 | * specifications of the PDB 3.2 format: http://www.wwpdb.org/docs.html
|
---|
| 382 | *
|
---|
| 383 | * A new atom is created and filled with available information, non-
|
---|
| 384 | * standard information is placed in additionalAtomData at the atom's id.
|
---|
[3ae731] | 385 | *
|
---|
| 386 | * \param line to parse as an atom
|
---|
[4fbca9c] | 387 | * \param newmol molecule to add parsed atoms to
|
---|
[3ae731] | 388 | */
|
---|
[4fbca9c] | 389 | void PdbParser::readAtomDataLine(std::string &line, molecule *newmol = NULL) {
|
---|
| 390 | vector<string>::iterator it;
|
---|
| 391 | stringstream lineStream;
|
---|
| 392 | atom* newAtom = World::getInstance().createAtom();
|
---|
[93fd43e] | 393 | PdbAtomInfoContainer &atomInfo = getadditionalAtomData(newAtom);
|
---|
[4fbca9c] | 394 | string word;
|
---|
| 395 | ConvertTo<size_t> toSize_t;
|
---|
| 396 | double tmp;
|
---|
| 397 |
|
---|
| 398 | lineStream << line;
|
---|
[16462f] | 399 | atomInfo.set(PdbKey::token, line.substr(0,6));
|
---|
[4fbca9c] | 400 | atomInfo.set(PdbKey::serial, line.substr(6,5));
|
---|
| 401 | std::pair< std::set<size_t>::const_iterator, bool> Inserter =
|
---|
[16462f] | 402 | SerialSet.insert(toSize_t(atomInfo.get<std::string>(PdbKey::serial)));
|
---|
[4fbca9c] | 403 | ASSERT(Inserter.second,
|
---|
| 404 | "PdbParser::readAtomDataLine() - ATOM contains entry with serial "
|
---|
[16462f] | 405 | +atomInfo.get<std::string>(PdbKey::serial)+" already present!");
|
---|
[4fbca9c] | 406 | // assign hightest+1 instead, but then beware of CONECT entries! Another map needed!
|
---|
| 407 | // if (!Inserter.second) {
|
---|
| 408 | // const size_t id = (*SerialSet.rbegin())+1;
|
---|
| 409 | // SerialSet.insert(id);
|
---|
| 410 | // atomInfo.set(PdbKey::serial, toString(id));
|
---|
| 411 | // DoeLog(2) && (eLog() << Verbose(2)
|
---|
[16462f] | 412 | // << "Serial " << atomInfo.get<std::string>(PdbKey::serial) << " already present, "
|
---|
[4fbca9c] | 413 | // << "assigning " << toString(id) << " instead." << std::endl);
|
---|
[bb6193] | 414 | // }
|
---|
[4fbca9c] | 415 |
|
---|
| 416 | // check whether serial exists, if so, assign next available
|
---|
| 417 |
|
---|
| 418 | // DoLog(2) && (Log() << Verbose(2) << "Split line:"
|
---|
| 419 | // << line.substr(6,5) << "|"
|
---|
| 420 | // << line.substr(12,4) << "|"
|
---|
| 421 | // << line.substr(16,1) << "|"
|
---|
| 422 | // << line.substr(17,3) << "|"
|
---|
| 423 | // << line.substr(21,1) << "|"
|
---|
| 424 | // << line.substr(22,4) << "|"
|
---|
| 425 | // << line.substr(26,1) << "|"
|
---|
| 426 | // << line.substr(30,8) << "|"
|
---|
| 427 | // << line.substr(38,8) << "|"
|
---|
| 428 | // << line.substr(46,8) << "|"
|
---|
| 429 | // << line.substr(54,6) << "|"
|
---|
| 430 | // << line.substr(60,6) << "|"
|
---|
| 431 | // << line.substr(76,2) << "|"
|
---|
| 432 | // << line.substr(78,2) << std::endl);
|
---|
| 433 |
|
---|
[16462f] | 434 | setSerial(toSize_t(atomInfo.get<std::string>(PdbKey::serial)), newAtom->getId());
|
---|
[4fbca9c] | 435 | atomInfo.set(PdbKey::name, line.substr(12,4));
|
---|
[16462f] | 436 | atomInfo.set(PdbKey::altLoc, line.substr(16,1));
|
---|
[4fbca9c] | 437 | atomInfo.set(PdbKey::resName, line.substr(17,3));
|
---|
| 438 | atomInfo.set(PdbKey::chainID, line.substr(21,1));
|
---|
| 439 | atomInfo.set(PdbKey::resSeq, line.substr(22,4));
|
---|
| 440 | atomInfo.set(PdbKey::iCode, line.substr(26,1));
|
---|
| 441 | PdbAtomInfoContainer::ScanKey(tmp, line.substr(30,8));
|
---|
| 442 | newAtom->set(0, tmp);
|
---|
| 443 | PdbAtomInfoContainer::ScanKey(tmp, line.substr(38,8));
|
---|
| 444 | newAtom->set(1, tmp);
|
---|
| 445 | PdbAtomInfoContainer::ScanKey(tmp, line.substr(46,8));
|
---|
| 446 | newAtom->set(2, tmp);
|
---|
| 447 | atomInfo.set(PdbKey::occupancy, line.substr(54,6));
|
---|
| 448 | atomInfo.set(PdbKey::tempFactor, line.substr(60,6));
|
---|
| 449 | atomInfo.set(PdbKey::charge, line.substr(78,2));
|
---|
[16462f] | 450 | atomInfo.set(PdbKey::element, line.substr(76,2));
|
---|
| 451 | const element *elem = World::getInstance().getPeriode()
|
---|
| 452 | ->FindElement(atomInfo.get<std::string>(PdbKey::element));
|
---|
| 453 | ASSERT(elem != NULL,
|
---|
| 454 | "PdbParser::readAtomDataLine() - element "+atomInfo.get<std::string>(PdbKey::element)+" is unknown!");
|
---|
| 455 | newAtom->setType(elem);
|
---|
[4fbca9c] | 456 |
|
---|
| 457 | if (newmol != NULL)
|
---|
| 458 | newmol->AddAtom(newAtom);
|
---|
| 459 |
|
---|
| 460 | // printAtomInfo(newAtom);
|
---|
[3ae731] | 461 | }
|
---|
| 462 |
|
---|
[4fbca9c] | 463 | /** Prints all PDB-specific information known about an atom.
|
---|
[3ae731] | 464 | *
|
---|
| 465 | */
|
---|
[4fbca9c] | 466 | void PdbParser::printAtomInfo(const atom * const newAtom) const
|
---|
| 467 | {
|
---|
| 468 | const PdbAtomInfoContainer &atomInfo = additionalAtomData.at(newAtom->getId()); // operator[] const does not exist
|
---|
| 469 |
|
---|
| 470 | DoLog(1) && (Log() << Verbose(1) << "We know about atom " << newAtom->getId() << ":" << std::endl);
|
---|
[16462f] | 471 | DoLog(1) && (Log() << Verbose(1) << "\ttoken is " << atomInfo.get<std::string>(PdbKey::token) << std::endl);
|
---|
| 472 | DoLog(1) && (Log() << Verbose(1) << "\tserial is " << atomInfo.get<int>(PdbKey::serial) << std::endl);
|
---|
| 473 | DoLog(1) && (Log() << Verbose(1) << "\tname is " << atomInfo.get<std::string>(PdbKey::name) << std::endl);
|
---|
| 474 | DoLog(1) && (Log() << Verbose(1) << "\taltLoc is " << atomInfo.get<std::string>(PdbKey::altLoc) << std::endl);
|
---|
| 475 | DoLog(1) && (Log() << Verbose(1) << "\tresName is " << atomInfo.get<std::string>(PdbKey::resName) << std::endl);
|
---|
| 476 | DoLog(1) && (Log() << Verbose(1) << "\tchainID is " << atomInfo.get<std::string>(PdbKey::chainID) << std::endl);
|
---|
| 477 | DoLog(1) && (Log() << Verbose(1) << "\tresSeq is " << atomInfo.get<int>(PdbKey::resSeq) << std::endl);
|
---|
| 478 | DoLog(1) && (Log() << Verbose(1) << "\tiCode is " << atomInfo.get<std::string>(PdbKey::iCode) << std::endl);
|
---|
| 479 | DoLog(1) && (Log() << Verbose(1) << "\tX is " << atomInfo.get<double>(PdbKey::X) << std::endl);
|
---|
| 480 | DoLog(1) && (Log() << Verbose(1) << "\tY is " << atomInfo.get<double>(PdbKey::Y) << std::endl);
|
---|
| 481 | DoLog(1) && (Log() << Verbose(1) << "\tZ is " << atomInfo.get<double>(PdbKey::Z) << std::endl);
|
---|
| 482 | DoLog(1) && (Log() << Verbose(1) << "\toccupancy is " << atomInfo.get<double>(PdbKey::occupancy) << std::endl);
|
---|
| 483 | DoLog(1) && (Log() << Verbose(1) << "\ttempFactor is " << atomInfo.get<double>(PdbKey::tempFactor) << std::endl);
|
---|
[4fbca9c] | 484 | DoLog(1) && (Log() << Verbose(1) << "\telement is '" << *(newAtom->getType()) << "'" << std::endl);
|
---|
[16462f] | 485 | DoLog(1) && (Log() << Verbose(1) << "\tcharge is " << atomInfo.get<int>(PdbKey::charge) << std::endl);
|
---|
[3ae731] | 486 | }
|
---|
| 487 |
|
---|
| 488 | /**
|
---|
[4fbca9c] | 489 | * Reads neighbor information for one atom from the input.
|
---|
| 490 | *
|
---|
| 491 | * \param line to parse as an atom
|
---|
[3ae731] | 492 | */
|
---|
[4fbca9c] | 493 | void PdbParser::readNeighbors(std::string &line)
|
---|
| 494 | {
|
---|
| 495 | const size_t length = line.length();
|
---|
| 496 | std::list<size_t> ListOfNeighbors;
|
---|
| 497 | ConvertTo<size_t> toSize_t;
|
---|
| 498 |
|
---|
| 499 | // obtain neighbours
|
---|
| 500 | // show split line for debugging
|
---|
| 501 | string output;
|
---|
| 502 | ASSERT(length >=16,
|
---|
| 503 | "PdbParser::readNeighbors() - CONECT entry has not enough entries: "+line+"!");
|
---|
| 504 | // output = "Split line:|";
|
---|
| 505 | // output += line.substr(6,5) + "|";
|
---|
| 506 | const size_t id = toSize_t(line.substr(6,5));
|
---|
| 507 | for (size_t index = 11; index <= 26; index+=5) {
|
---|
| 508 | if (index+5 <= length) {
|
---|
| 509 | // output += line.substr(index,5) + "|";
|
---|
| 510 | const size_t otherid = toSize_t(line.substr(index,5));
|
---|
| 511 | ListOfNeighbors.push_back(otherid);
|
---|
| 512 | } else {
|
---|
| 513 | break;
|
---|
| 514 | }
|
---|
| 515 | }
|
---|
| 516 | // DoLog(2) && (Log() << Verbose(2) << output << std::endl);
|
---|
| 517 |
|
---|
| 518 | // add neighbours
|
---|
[16462f] | 519 | atom *_atom = World::getInstance().getAtom(AtomById(getSerial(id)));
|
---|
[4fbca9c] | 520 | for (std::list<size_t>::const_iterator iter = ListOfNeighbors.begin();
|
---|
| 521 | iter != ListOfNeighbors.end();
|
---|
| 522 | ++iter) {
|
---|
| 523 | // DoLog(1) && (Log() << Verbose(1) << "Adding Bond (" << getAtomId(id) << "," << getAtomId(*iter) << ")" << std::endl);
|
---|
[16462f] | 524 | atom * const _Otheratom = World::getInstance().getAtom(AtomById(getSerial(*iter)));
|
---|
[073a9e4] | 525 | _atom->addBond(WorldTime::getTime(), _Otheratom);
|
---|
[4fbca9c] | 526 | }
|
---|
[3ae731] | 527 | }
|
---|
| 528 |
|
---|
| 529 | /**
|
---|
| 530 | * Replaces atom IDs read from the file by the corresponding world IDs. All IDs
|
---|
| 531 | * IDs of the input string will be replaced; expected separating characters are
|
---|
| 532 | * "-" and ",".
|
---|
| 533 | *
|
---|
| 534 | * \param string in which atom IDs should be adapted
|
---|
| 535 | *
|
---|
| 536 | * \return input string with modified atom IDs
|
---|
| 537 | */
|
---|
[4fbca9c] | 538 | //string PdbParser::adaptIdDependentDataString(string data) {
|
---|
[bb6193] | 539 | // // there might be no IDs
|
---|
| 540 | // if (data == "-") {
|
---|
| 541 | // return "-";
|
---|
| 542 | // }
|
---|
| 543 | //
|
---|
| 544 | // char separator;
|
---|
| 545 | // int id;
|
---|
| 546 | // stringstream line, result;
|
---|
| 547 | // line << data;
|
---|
| 548 | //
|
---|
| 549 | // line >> id;
|
---|
| 550 | // result << atomIdMap[id];
|
---|
| 551 | // while (line.good()) {
|
---|
| 552 | // line >> separator >> id;
|
---|
| 553 | // result << separator << atomIdMap[id];
|
---|
| 554 | // }
|
---|
| 555 | //
|
---|
| 556 | // return result.str();
|
---|
[4fbca9c] | 557 | // return "";
|
---|
| 558 | //}
|
---|
[3ae731] | 559 |
|
---|
| 560 |
|
---|
[4fbca9c] | 561 | bool PdbParser::operator==(const PdbParser& b) const
|
---|
| 562 | {
|
---|
| 563 | bool status = true;
|
---|
| 564 | World::AtomComposite atoms = World::getInstance().getAllAtoms();
|
---|
| 565 | for (World::AtomComposite::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
|
---|
| 566 | if ((additionalAtomData.find((*iter)->getId()) != additionalAtomData.end())
|
---|
| 567 | && (b.additionalAtomData.find((*iter)->getId()) != b.additionalAtomData.end())) {
|
---|
| 568 | const PdbAtomInfoContainer &atomInfo = additionalAtomData.at((*iter)->getId());
|
---|
| 569 | const PdbAtomInfoContainer &OtheratomInfo = b.additionalAtomData.at((*iter)->getId());
|
---|
| 570 |
|
---|
[16462f] | 571 | status = status && (atomInfo.get<std::string>(PdbKey::serial) == OtheratomInfo.get<std::string>(PdbKey::serial));
|
---|
[4fbca9c] | 572 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in serials!" << std::endl);
|
---|
[16462f] | 573 | status = status && (atomInfo.get<std::string>(PdbKey::name) == OtheratomInfo.get<std::string>(PdbKey::name));
|
---|
[4fbca9c] | 574 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in names!" << std::endl);
|
---|
[16462f] | 575 | status = status && (atomInfo.get<std::string>(PdbKey::altLoc) == OtheratomInfo.get<std::string>(PdbKey::altLoc));
|
---|
| 576 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in altLocs!" << std::endl);
|
---|
| 577 | status = status && (atomInfo.get<std::string>(PdbKey::resName) == OtheratomInfo.get<std::string>(PdbKey::resName));
|
---|
[4fbca9c] | 578 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resNames!" << std::endl);
|
---|
[16462f] | 579 | status = status && (atomInfo.get<std::string>(PdbKey::chainID) == OtheratomInfo.get<std::string>(PdbKey::chainID));
|
---|
[4fbca9c] | 580 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in chainIDs!" << std::endl);
|
---|
[16462f] | 581 | status = status && (atomInfo.get<std::string>(PdbKey::resSeq) == OtheratomInfo.get<std::string>(PdbKey::resSeq));
|
---|
[4fbca9c] | 582 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resSeqs!" << std::endl);
|
---|
[16462f] | 583 | status = status && (atomInfo.get<std::string>(PdbKey::iCode) == OtheratomInfo.get<std::string>(PdbKey::iCode));
|
---|
[4fbca9c] | 584 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in iCodes!" << std::endl);
|
---|
[16462f] | 585 | status = status && (atomInfo.get<std::string>(PdbKey::occupancy) == OtheratomInfo.get<std::string>(PdbKey::occupancy));
|
---|
[4fbca9c] | 586 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in occupancies!" << std::endl);
|
---|
[16462f] | 587 | status = status && (atomInfo.get<std::string>(PdbKey::tempFactor) == OtheratomInfo.get<std::string>(PdbKey::tempFactor));
|
---|
[4fbca9c] | 588 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in tempFactors!" << std::endl);
|
---|
[16462f] | 589 | status = status && (atomInfo.get<std::string>(PdbKey::charge) == OtheratomInfo.get<std::string>(PdbKey::charge));
|
---|
[4fbca9c] | 590 | if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in charges!" << std::endl);
|
---|
| 591 | }
|
---|
[3ae731] | 592 | }
|
---|
| 593 |
|
---|
[4fbca9c] | 594 | return status;
|
---|
[3ae731] | 595 | }
|
---|
| 596 |
|
---|