source: src/Parser/PdbParser.cpp@ bb6193

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since bb6193 was bb6193, checked in by Frederik Heber <heber@…>, 14 years ago

Added PdbParser with save capability.

  • load does not yet work.
  • added test part to Simple_configuration/2 (new file test.pdb)
  • Property mode set to 100644
File size: 11.8 KB
Line 
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
20#include "Helpers/MemDebug.hpp"
21
22#include "Helpers/Assert.hpp"
23#include "Helpers/Log.hpp"
24#include "Helpers/Verbose.hpp"
25#include "PdbParser.hpp"
26#include "World.hpp"
27#include "atom.hpp"
28#include "bond.hpp"
29#include "element.hpp"
30#include "molecule.hpp"
31#include "periodentafel.hpp"
32#include "Descriptors/AtomIdDescriptor.hpp"
33
34#include <map>
35#include <vector>
36
37#include <iostream>
38#include <iomanip>
39
40using namespace std;
41
42/**
43 * Constructor.
44 */
45PdbParser::PdbParser() {
46 knownKeys[" "] = PdbKey::noKey; // with this we can detect invalid keys
47 knownKeys["x"] = PdbKey::x;
48 knownKeys["Id"] = PdbKey::Id;
49 knownKeys["Type"] = PdbKey::Type;
50 knownKeys["extType"] = PdbKey::extType;
51 knownKeys["name"] = PdbKey::name;
52 knownKeys["resName"] = PdbKey::resName;
53 knownKeys["chainID"] = PdbKey::chainID;
54 knownKeys["resSeq"] = PdbKey::resSeq;
55 knownKeys["occupancy"] = PdbKey::occupancy;
56 knownKeys["tempFactor"] = PdbKey::tempFactor;
57 knownKeys["segID"] = PdbKey::segID;
58 knownKeys["charge"] = PdbKey::charge;
59}
60
61/**
62 * Destructor.
63 */
64PdbParser::~PdbParser() {
65 additionalAtomData.clear();
66 atomIdMap.clear();
67 knownKeys.clear();
68}
69
70/**
71 * Loads atoms from a tremolo-formatted file.
72 *
73 * \param tremolo file
74 */
75void PdbParser::load(istream* file) {
76// string line;
77// string::size_type location;
78//
79// usedFields.clear();
80// while (file->good()) {
81// std::getline(*file, line, '\n');
82// if (usedFields.empty()) {
83// location = line.find("ATOMDATA", 0);
84// if (location != string::npos) {
85// parseAtomDataKeysLine(line, location + 8);
86// }
87// }
88// if (line.length() > 0 && line.at(0) != '#') {
89// readAtomDataLine(line);
90// }
91// }
92//
93// processNeighborInformation();
94// adaptImprData();
95// adaptTorsion();
96}
97
98/**
99 * Saves the World's current state into as a tremolo file.
100 *
101 * \param file where to save the state
102 */
103void PdbParser::save(ostream* file) {
104 DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl);
105
106 {
107 // add initial remark
108 *file << "REMARK created by molecuilder on ";
109 time_t now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time'
110 // ctime ends in \n\0, we have to cut away the newline
111 std::string time(ctime(&now));
112 size_t pos = time.find('\n');
113 if (pos != 0)
114 *file << time.substr(0,pos);
115 else
116 *file << time;
117 *file << endl;
118 }
119
120 {
121 vector<atom *> AtomList = World::getInstance().getAllAtoms();
122
123 std::vector<int> elementNo(MAX_ELEMENTS,1);
124 char name[MAXSTRINGSIZE];
125
126 // write ATOMs
127 int AtomNo = 1; // serial number starts at 1 in pdb
128 int MolNo = 1; // residue number starts at 1 in pdb
129 for (vector<atom *>::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
130 const size_t Z = (*atomIt)->getType()->getAtomicNumber();
131 sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), elementNo[Z]);
132 elementNo[Z] = (elementNo[Z]+1) % 100; // confine to two digits
133 const molecule *mol = (*atomIt)->getMolecule();
134 if (mol == NULL) { // for homeless atoms, MolNo = 0 is reserved
135 MolNo = 0;
136 } else {
137 MolNo = mol->getId();
138 }
139 saveLine(file, *atomIt, name, AtomNo, MolNo);
140 atomIdMap.insert( pair<int, int>((*atomIt)->getId(), AtomNo) );
141 AtomNo++;
142 }
143
144 // write CONECTs
145 for (vector<atom *>::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
146 writeNeighbors(file, 4, *atomIt);
147 }
148 }
149
150 // END
151 *file << "END" << endl;
152}
153
154/**
155 * Writes one line of tremolo-formatted data to the provided stream.
156 *
157 * \param stream where to write the line to
158 * \param *currentAtom the atom of which information should be written
159 * \param *name name of atom, i.e. H01
160 * \param AtomNo serial number of atom
161 * \param ResidueNo number of residue
162 */
163void PdbParser::saveLine(ostream* file, const atom* currentAtom, const char *name, const int AtomNo, const int ResidueNo) {
164 *file << "ATOM ";
165 *file << setw(6) << AtomNo; /* atom serial number */
166 *file << setw(1) << " ";
167 *file << setfill(' ') << left << setw(4) << name << right; /* atom name */
168 *file << setw(1) << " ";
169 *file << setfill(' ') << setw(3) << ((currentAtom->getMolecule() != NULL) ? currentAtom->getMolecule()->getName().substr(0,3) : "-"); /* residue name */
170 *file << setw(1) << " ";
171 *file << setfill(' ') << setw(1) << (char)('a'+(unsigned char)(AtomNo % 26)); /* letter for chain */
172 *file << setw(4) << ResidueNo; /* residue sequence number */
173 *file << setw(4) << " ";
174 for (int i=0;i<NDIM;i++) {
175 *file << setw(8) << setprecision(3) << showpoint << currentAtom->at(i); /* positional coordinate in Angstroem */
176 }
177 *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getValence(); /* occupancy */
178 *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getNoValenceOrbitals(); /* temperature factor */
179 *file << noshowpoint;
180 *file << setw(6) << " ";
181 *file << setw(4) << "0";
182 *file << setfill(' ') << setw(2) << currentAtom->getType()->getSymbol();
183 *file << setw(2) << "0";
184
185 *file << endl;
186}
187
188/**
189 * Writes the neighbor information of one atom to the provided stream.
190 *
191 * \param *file where to write neighbor information to
192 * \param MaxnumberOfNeighbors of neighbors
193 * \param *currentAtom to the atom of which to take the neighbor information
194 */
195void PdbParser::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) {
196 if (!currentAtom->ListOfBonds.empty()) {
197 *file << "CONECT";
198 int MaxNo = 0;
199 for(BondList::iterator currentBond = currentAtom->ListOfBonds.begin(); currentBond != currentAtom->ListOfBonds.end(); ++currentBond) {
200 if (MaxNo < MaxnumberOfNeighbors) {
201 ASSERT(atomIdMap.find((*currentBond)->GetOtherAtom(currentAtom)->getId()) != atomIdMap.end(), "Id of atom not stored in PdbParser::atomIdMap.");
202 *file << setw(5) << atomIdMap[(*currentBond)->GetOtherAtom(currentAtom)->getId()];
203 }
204 MaxNo++;
205 }
206 *file << endl;
207 }
208}
209
210/**
211 * Reads one data line of a tremolo file and interprets it according to the keys
212 * obtained from the ATOMDATA line.
213 *
214 * \param line to parse as an atom
215 */
216void PdbParser::readAtomDataLine(string line) {
217// vector<string>::iterator it;
218// stringstream lineStream;
219// atom* newAtom = World::getInstance().createAtom();
220// PdbAtomInfoContainer *atomInfo = NULL;
221// additionalAtomData[newAtom->getId()] = *(new PdbAtomInfoContainer);
222// atomInfo = &additionalAtomData[newAtom->getId()];
223// PdbKey::atomDataKey currentField;
224// string word;
225// int oldId;
226// double tmp;
227//
228// lineStream << line;
229// for (it = usedFields.begin(); it < usedFields.end(); it++) {
230// currentField = knownKeys[it->substr(0, it->find("="))];
231// switch (currentField) {
232// case PdbKey::x :
233// // for the moment, assume there are always three dimensions
234// for (int i=0;i<NDIM;i++) {
235// lineStream >> tmp;
236// newAtom->set(i, tmp);
237// }
238// break;
239// case PdbKey::u :
240// // for the moment, assume there are always three dimensions
241// lineStream >> newAtom->AtomicVelocity[0];
242// lineStream >> newAtom->AtomicVelocity[1];
243// lineStream >> newAtom->AtomicVelocity[2];
244// break;
245// case PdbKey::Type :
246// char type[3];
247// lineStream >> type;
248// newAtom->setType(World::getInstance().getPeriode()->FindElement(type));
249// ASSERT(newAtom->getType(), "Type was not set for this atom");
250// break;
251// case PdbKey::Id :
252// lineStream >> oldId;
253// atomIdMap[oldId] = newAtom->getId();
254// break;
255// case PdbKey::neighbors :
256// readNeighbors(&lineStream,
257// atoi(it->substr(it->find("=") + 1, 1).c_str()), newAtom->getId());
258// break;
259// default :
260// lineStream >> word;
261// atomInfo->set(currentField, word);
262// break;
263// }
264// }
265}
266
267/**
268 * Reads neighbor information for one atom from the input.
269 *
270 * \param stream where to read the information from
271 * \param number of neighbors to read
272 * \param world id of the atom the information belongs to
273 */
274void PdbParser::readNeighbors(stringstream* line, int numberOfNeighbors, int atomId) {
275// int neighborId = 0;
276// for (int i = 0; i < numberOfNeighbors; i++) {
277// *line >> neighborId;
278// // 0 is used to fill empty neighbor positions in the tremolo file.
279// if (neighborId > 0) {
280// additionalAtomData[atomId].neighbors.push_back(neighborId);
281// }
282// }
283}
284
285/**
286 * Adds the collected neighbor information to the atoms in the world. The atoms
287 * are found by their current ID and mapped to the corresponding atoms with the
288 * Id found in the parsed file.
289 */
290void PdbParser::processNeighborInformation() {
291// if (!isUsedField("neighbors")) {
292// return;
293// }
294//
295// for(map<int, PdbAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
296// currentInfo != additionalAtomData.end(); currentInfo++
297// ) {
298// for(vector<int>::iterator neighbor = currentInfo->second.neighbors.begin();
299// neighbor != currentInfo->second.neighbors.end(); neighbor++
300// ) {
301// World::getInstance().getAtom(AtomById(currentInfo->first))
302// ->addBond(World::getInstance().getAtom(AtomById(atomIdMap[*neighbor])));
303// }
304// }
305}
306
307/**
308 * Replaces atom IDs read from the file by the corresponding world IDs. All IDs
309 * IDs of the input string will be replaced; expected separating characters are
310 * "-" and ",".
311 *
312 * \param string in which atom IDs should be adapted
313 *
314 * \return input string with modified atom IDs
315 */
316string PdbParser::adaptIdDependentDataString(string data) {
317// // there might be no IDs
318// if (data == "-") {
319// return "-";
320// }
321//
322// char separator;
323// int id;
324// stringstream line, result;
325// line << data;
326//
327// line >> id;
328// result << atomIdMap[id];
329// while (line.good()) {
330// line >> separator >> id;
331// result << separator << atomIdMap[id];
332// }
333//
334// return result.str();
335 return "";
336}
337
338
339PdbAtomInfoContainer::PdbAtomInfoContainer() :
340 name("-"),
341 resName("-"),
342 chainID("0"),
343 resSeq("0"),
344 occupancy("0"),
345 tempFactor("0"),
346 segID("0"),
347 charge("0")
348{}
349
350void PdbAtomInfoContainer::set(PdbKey::PdbDataKey key, string value) {
351 switch (key) {
352 case PdbKey::extType :
353 extType = value;
354 break;
355 case PdbKey::name :
356 name = value;
357 break;
358 case PdbKey::resName :
359 resName = value;
360 break;
361 case PdbKey::chainID :
362 chainID = value;
363 break;
364 case PdbKey::resSeq :
365 resSeq = value;
366 break;
367 case PdbKey::occupancy :
368 occupancy = value;
369 break;
370 case PdbKey::tempFactor :
371 tempFactor = value;
372 break;
373 case PdbKey::segID :
374 segID = value;
375 break;
376 case PdbKey::charge :
377 charge = value;
378 break;
379 default :
380 cout << "Unknown key: " << key << ", value: " << value << endl;
381 break;
382 }
383}
384
385string PdbAtomInfoContainer::get(PdbKey::PdbDataKey key) {
386 switch (key) {
387 case PdbKey::extType :
388 return extType;
389 case PdbKey::name :
390 return name;
391 case PdbKey::resName :
392 return resName;
393 case PdbKey::chainID :
394 return chainID;
395 case PdbKey::resSeq :
396 return resSeq;
397 case PdbKey::occupancy :
398 return occupancy;
399 case PdbKey::tempFactor :
400 return tempFactor;
401 case PdbKey::segID :
402 return segID;
403 case PdbKey::charge :
404 return charge;
405 default :
406 cout << "Unknown key: " << key << endl;
407 return "";
408 }
409}
410
Note: See TracBrowser for help on using the repository browser.