| [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 |  | 
|---|