1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
|
---|
5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
6 | */
|
---|
7 |
|
---|
8 | /*
|
---|
9 | * TremoloParser.cpp
|
---|
10 | *
|
---|
11 | * Created on: Mar 2, 2010
|
---|
12 | * Author: metzler
|
---|
13 | */
|
---|
14 |
|
---|
15 | // include config.h
|
---|
16 | #ifdef HAVE_CONFIG_H
|
---|
17 | #include <config.h>
|
---|
18 | #endif
|
---|
19 |
|
---|
20 | #include "CodePatterns/MemDebug.hpp"
|
---|
21 |
|
---|
22 | #include "CodePatterns/Assert.hpp"
|
---|
23 | #include "CodePatterns/Log.hpp"
|
---|
24 | #include "CodePatterns/toString.hpp"
|
---|
25 | #include "CodePatterns/Verbose.hpp"
|
---|
26 |
|
---|
27 | #include "TremoloParser.hpp"
|
---|
28 |
|
---|
29 | #include "Atom/atom.hpp"
|
---|
30 | #include "Bond/bond.hpp"
|
---|
31 | #include "Box.hpp"
|
---|
32 | #include "Descriptors/AtomIdDescriptor.hpp"
|
---|
33 | #include "Element/element.hpp"
|
---|
34 | #include "Element/periodentafel.hpp"
|
---|
35 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
36 | #include "molecule.hpp"
|
---|
37 | #include "MoleculeListClass.hpp"
|
---|
38 | #include "World.hpp"
|
---|
39 | #include "WorldTime.hpp"
|
---|
40 |
|
---|
41 |
|
---|
42 | #include <algorithm>
|
---|
43 | #include <boost/lexical_cast.hpp>
|
---|
44 | #include <boost/tokenizer.hpp>
|
---|
45 | #include <iostream>
|
---|
46 | #include <iomanip>
|
---|
47 | #include <map>
|
---|
48 | #include <sstream>
|
---|
49 | #include <vector>
|
---|
50 |
|
---|
51 | // declare specialized static variables
|
---|
52 | const std::string FormatParserTrait<tremolo>::name = "tremolo";
|
---|
53 | const std::string FormatParserTrait<tremolo>::suffix = "data";
|
---|
54 | const ParserTypes FormatParserTrait<tremolo>::type = tremolo;
|
---|
55 |
|
---|
56 | /**
|
---|
57 | * Constructor.
|
---|
58 | */
|
---|
59 | FormatParser< tremolo >::FormatParser() :
|
---|
60 | FormatParser_common(NULL)
|
---|
61 | {
|
---|
62 | knownKeys["x"] = TremoloKey::x;
|
---|
63 | knownKeys["u"] = TremoloKey::u;
|
---|
64 | knownKeys["F"] = TremoloKey::F;
|
---|
65 | knownKeys["stress"] = TremoloKey::stress;
|
---|
66 | knownKeys["Id"] = TremoloKey::Id;
|
---|
67 | knownKeys["neighbors"] = TremoloKey::neighbors;
|
---|
68 | knownKeys["imprData"] = TremoloKey::imprData;
|
---|
69 | knownKeys["GroupMeasureTypeNo"] = TremoloKey::GroupMeasureTypeNo;
|
---|
70 | knownKeys["type"] = TremoloKey::type;
|
---|
71 | knownKeys["extType"] = TremoloKey::extType;
|
---|
72 | knownKeys["name"] = TremoloKey::name;
|
---|
73 | knownKeys["resName"] = TremoloKey::resName;
|
---|
74 | knownKeys["chainID"] = TremoloKey::chainID;
|
---|
75 | knownKeys["resSeq"] = TremoloKey::resSeq;
|
---|
76 | knownKeys["occupancy"] = TremoloKey::occupancy;
|
---|
77 | knownKeys["tempFactor"] = TremoloKey::tempFactor;
|
---|
78 | knownKeys["segID"] = TremoloKey::segID;
|
---|
79 | knownKeys["Charge"] = TremoloKey::Charge;
|
---|
80 | knownKeys["charge"] = TremoloKey::charge;
|
---|
81 | knownKeys["GrpTypeNo"] = TremoloKey::GrpTypeNo;
|
---|
82 | knownKeys["torsion"] = TremoloKey::torsion;
|
---|
83 | knownKeys[" "] = TremoloKey::noKey; // with this we can detect invalid keys
|
---|
84 |
|
---|
85 | createKnownTypesByIdentity();
|
---|
86 |
|
---|
87 | // invert knownKeys for debug output
|
---|
88 | for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin(); iter != knownKeys.end(); ++iter)
|
---|
89 | knownKeyNames.insert( make_pair( iter->second, iter->first) );
|
---|
90 |
|
---|
91 | additionalAtomData.clear();
|
---|
92 | }
|
---|
93 |
|
---|
94 | /**
|
---|
95 | * Destructor.
|
---|
96 | */
|
---|
97 | FormatParser< tremolo >::~FormatParser()
|
---|
98 | {
|
---|
99 | usedFields_save.clear();
|
---|
100 | additionalAtomData.clear();
|
---|
101 | knownKeys.clear();
|
---|
102 | }
|
---|
103 |
|
---|
104 | /**
|
---|
105 | * Loads atoms from a tremolo-formatted file.
|
---|
106 | *
|
---|
107 | * \param tremolo file
|
---|
108 | */
|
---|
109 | void FormatParser< tremolo >::load(istream* file) {
|
---|
110 | std::string line;
|
---|
111 | std::string::size_type location;
|
---|
112 |
|
---|
113 | // reset the id maps
|
---|
114 | resetIdAssociations();
|
---|
115 |
|
---|
116 | molecule *newmol = World::getInstance().createMolecule();
|
---|
117 | newmol->ActiveFlag = true;
|
---|
118 | // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
|
---|
119 | World::getInstance().getMolecules()->insert(newmol);
|
---|
120 | while (file->good()) {
|
---|
121 | std::getline(*file, line, '\n');
|
---|
122 | // we only parse in the first ATOMDATA line
|
---|
123 | if (usedFields_load.empty()) {
|
---|
124 | location = line.find("ATOMDATA", 0);
|
---|
125 | if (location != string::npos) {
|
---|
126 | parseAtomDataKeysLine(line, location + 8, usedFields_load);
|
---|
127 | }
|
---|
128 | }
|
---|
129 | if (line.length() > 0 && line.at(0) != '#') {
|
---|
130 | readAtomDataLine(line, newmol);
|
---|
131 | }
|
---|
132 | }
|
---|
133 | LOG(3, "DEBUG: Local usedFields is: " << usedFields_load);
|
---|
134 |
|
---|
135 | // refresh atom::nr and atom::name
|
---|
136 | std::vector<atomId_t> atoms(newmol->getAtomCount());
|
---|
137 | std::transform(newmol->begin(), newmol->end(), atoms.begin(),
|
---|
138 | boost::bind(&atom::getId, _1));
|
---|
139 | processNeighborInformation(atoms);
|
---|
140 | adaptImprData();
|
---|
141 | adaptTorsion();
|
---|
142 |
|
---|
143 | // append usedFields to global usedFields, is made unique on save, clear after use
|
---|
144 | usedFields_save.insert(usedFields_save.end(), usedFields_load.begin(), usedFields_load.end());
|
---|
145 | usedFields_load.clear();
|
---|
146 | }
|
---|
147 |
|
---|
148 | /**
|
---|
149 | * Saves the \a atoms into as a tremolo file.
|
---|
150 | *
|
---|
151 | * \param file where to save the state
|
---|
152 | * \param atoms atoms to store
|
---|
153 | */
|
---|
154 | void FormatParser< tremolo >::save(std::ostream* file, const std::vector<atom *> &AtomList) {
|
---|
155 | LOG(0, "Saving changes to tremolo.");
|
---|
156 |
|
---|
157 | // install default usedFields if empty so far
|
---|
158 | if (usedFields_save.empty()) {
|
---|
159 | // default behavior: use all possible keys on output
|
---|
160 | for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin();
|
---|
161 | iter != knownKeys.end(); ++iter)
|
---|
162 | if (iter->second != TremoloKey::noKey) // don't add noKey
|
---|
163 | usedFields_save.push_back(iter->first);
|
---|
164 | }
|
---|
165 | // make present usedFields_save unique
|
---|
166 | makeUsedFieldsUnique(usedFields_save);
|
---|
167 | LOG(1, "DEBUG: Global (with unique entries) usedFields_save is: " << usedFields_save);
|
---|
168 |
|
---|
169 | // distribute ids continuously
|
---|
170 | distributeContinuousIds(AtomList);
|
---|
171 |
|
---|
172 | // store atomdata
|
---|
173 | save_AtomDataLine(file);
|
---|
174 |
|
---|
175 | // store box
|
---|
176 | save_BoxLine(file);
|
---|
177 |
|
---|
178 | // store particles
|
---|
179 | for (std::vector<atom*>::const_iterator atomIt = AtomList.begin();
|
---|
180 | atomIt != AtomList.end(); ++atomIt)
|
---|
181 | saveLine(file, *atomIt);
|
---|
182 | }
|
---|
183 |
|
---|
184 | /** Helper function to make \given fields unique while preserving the order of first appearance.
|
---|
185 | *
|
---|
186 | * As std::unique only removes element if equal to predecessor, a vector is only
|
---|
187 | * made unique if sorted beforehand. But sorting would destroy order of first
|
---|
188 | * appearance, hence we do the sorting on a temporary field and add the unique
|
---|
189 | * elements in the order as in \a fields.
|
---|
190 | *
|
---|
191 | * @param fields usedFields to make unique while preserving order of appearance
|
---|
192 | */
|
---|
193 | void FormatParser< tremolo >::makeUsedFieldsUnique(usedFields_t &fields)
|
---|
194 | {
|
---|
195 | // std::unique only removes if predecessor is equal, not over whole range, hence do it manually
|
---|
196 | usedFields_t temp_fields(usedFields_save);
|
---|
197 | std::sort(temp_fields.begin(), temp_fields.end());
|
---|
198 | usedFields_t::iterator it =
|
---|
199 | std::unique(temp_fields.begin(), temp_fields.end()); // skips all duplicates in the vector
|
---|
200 | temp_fields.erase(it, temp_fields.end());
|
---|
201 | usedFields_t usedfields(usedFields_save);
|
---|
202 | usedFields_save.clear();
|
---|
203 | usedFields_save.reserve(temp_fields.size());
|
---|
204 | // now go through each usedFields entry, check if in temp_fields and remove there on first occurence
|
---|
205 | for (usedFields_t::const_iterator iter = usedfields.begin();
|
---|
206 | iter != usedfields.end(); ++iter) {
|
---|
207 | usedFields_t::iterator uniqueiter =
|
---|
208 | std::find(temp_fields.begin(), temp_fields.end(), *iter);
|
---|
209 | if (uniqueiter != temp_fields.end()) {
|
---|
210 | usedFields_save.push_back(*iter);
|
---|
211 | // add only once to ATOMDATA
|
---|
212 | temp_fields.erase(uniqueiter);
|
---|
213 | }
|
---|
214 | }
|
---|
215 | ASSERT( temp_fields.empty(),
|
---|
216 | "FormatParser< tremolo >::save() - still unique entries left in temp_fields after unique?");
|
---|
217 | }
|
---|
218 |
|
---|
219 |
|
---|
220 | /** Resets and distributes the indices continuously.
|
---|
221 | *
|
---|
222 | * \param atoms atoms to store
|
---|
223 | */
|
---|
224 | void FormatParser< tremolo >::distributeContinuousIds(const std::vector<atom *> &AtomList)
|
---|
225 | {
|
---|
226 | resetIdAssociations();
|
---|
227 | atomId_t lastid = 0;
|
---|
228 | for (std::vector<atom*>::const_iterator atomIt = AtomList.begin();
|
---|
229 | atomIt != AtomList.end(); ++atomIt)
|
---|
230 | associateLocaltoGlobalId(++lastid, (*atomIt)->getId());
|
---|
231 | }
|
---|
232 |
|
---|
233 | /** Store Atomdata line to \a file.
|
---|
234 | *
|
---|
235 | * @param file output stream
|
---|
236 | */
|
---|
237 | void FormatParser< tremolo >::save_AtomDataLine(std::ostream* file) const
|
---|
238 | {
|
---|
239 | *file << "# ATOMDATA";
|
---|
240 | for (usedFields_t::const_iterator it=usedFields_save.begin();
|
---|
241 | it != usedFields_save.end(); ++it)
|
---|
242 | *file << "\t" << *it;
|
---|
243 | *file << std::endl;
|
---|
244 | }
|
---|
245 |
|
---|
246 | /** Store Box info to \a file
|
---|
247 | *
|
---|
248 | * @param file output stream
|
---|
249 | */
|
---|
250 | void FormatParser< tremolo >::save_BoxLine(std::ostream* file) const
|
---|
251 | {
|
---|
252 | *file << "# Box";
|
---|
253 | const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
|
---|
254 | for (size_t i=0; i<NDIM;++i)
|
---|
255 | for (size_t j=0; j<NDIM;++j)
|
---|
256 | *file << "\t" << M.at(i,j);
|
---|
257 | *file << std::endl;
|
---|
258 | }
|
---|
259 |
|
---|
260 | /** Add default info, when new atom is added to World.
|
---|
261 | *
|
---|
262 | * @param id of atom
|
---|
263 | */
|
---|
264 | void FormatParser< tremolo >::AtomInserted(atomId_t id)
|
---|
265 | {
|
---|
266 | std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
|
---|
267 | ASSERT(iter == additionalAtomData.end(),
|
---|
268 | "FormatParser< tremolo >::AtomInserted() - additionalAtomData already present for newly added atom "
|
---|
269 | +toString(id)+".");
|
---|
270 | // don't add entry, as this gives a default resSeq of 0 not the molecule id
|
---|
271 | // additionalAtomData.insert( std::make_pair(id, TremoloAtomInfoContainer()) );
|
---|
272 | }
|
---|
273 |
|
---|
274 | /** Remove additional AtomData info, when atom has been removed from World.
|
---|
275 | *
|
---|
276 | * @param id of atom
|
---|
277 | */
|
---|
278 | void FormatParser< tremolo >::AtomRemoved(atomId_t id)
|
---|
279 | {
|
---|
280 | std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
|
---|
281 | // as we do not insert AtomData on AtomInserted, we cannot be assured of its presence
|
---|
282 | // ASSERT(iter != additionalAtomData.end(),
|
---|
283 | // "FormatParser< tremolo >::AtomRemoved() - additionalAtomData is not present for atom "
|
---|
284 | // +toString(id)+" to remove.");
|
---|
285 | if (iter != additionalAtomData.end())
|
---|
286 | additionalAtomData.erase(iter);
|
---|
287 | }
|
---|
288 |
|
---|
289 | /**
|
---|
290 | * Writes one line of tremolo-formatted data to the provided stream.
|
---|
291 | *
|
---|
292 | * \param stream where to write the line to
|
---|
293 | * \param reference to the atom of which information should be written
|
---|
294 | */
|
---|
295 | void FormatParser< tremolo >::saveLine(std::ostream* file, const atom* currentAtom)
|
---|
296 | {
|
---|
297 | TremoloKey::atomDataKey currentField;
|
---|
298 |
|
---|
299 | LOG(4, "INFO: Saving atom " << *currentAtom << ", its father id is " << currentAtom->GetTrueFather()->getId());
|
---|
300 |
|
---|
301 | for (usedFields_t::iterator it = usedFields_save.begin(); it != usedFields_save.end(); it++) {
|
---|
302 | currentField = knownKeys[it->substr(0, it->find("="))];
|
---|
303 | switch (currentField) {
|
---|
304 | case TremoloKey::x :
|
---|
305 | // for the moment, assume there are always three dimensions
|
---|
306 | LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getPosition());
|
---|
307 | *file << currentAtom->at(0) << "\t";
|
---|
308 | *file << currentAtom->at(1) << "\t";
|
---|
309 | *file << currentAtom->at(2) << "\t";
|
---|
310 | break;
|
---|
311 | case TremoloKey::u :
|
---|
312 | // for the moment, assume there are always three dimensions
|
---|
313 | LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicVelocity());
|
---|
314 | *file << currentAtom->getAtomicVelocity()[0] << "\t";
|
---|
315 | *file << currentAtom->getAtomicVelocity()[1] << "\t";
|
---|
316 | *file << currentAtom->getAtomicVelocity()[2] << "\t";
|
---|
317 | break;
|
---|
318 | case TremoloKey::type :
|
---|
319 | if (additionalAtomData.count(currentAtom->getId())) {
|
---|
320 | if (additionalAtomData[currentAtom->getId()].get(currentField) != "-") {
|
---|
321 | LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
|
---|
322 | *file << additionalAtomData[currentAtom->getId()].get(currentField) << "\t";
|
---|
323 | } else {
|
---|
324 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value: " << currentAtom->getType()->getSymbol());
|
---|
325 | *file << currentAtom->getType()->getSymbol() << "\t";
|
---|
326 | }
|
---|
327 | } else if (additionalAtomData.count(currentAtom->GetTrueFather()->getId())) {
|
---|
328 | if (additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField) != "-") {
|
---|
329 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " stuff from father: " << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField));
|
---|
330 | *file << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField) << "\t";
|
---|
331 | } else {
|
---|
332 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value from father: " << currentAtom->GetTrueFather()->getType()->getSymbol());
|
---|
333 | *file << currentAtom->GetTrueFather()->getType()->getSymbol() << "\t";
|
---|
334 | }
|
---|
335 | } else {
|
---|
336 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " its default value: " << currentAtom->getType()->getSymbol());
|
---|
337 | *file << currentAtom->getType()->getSymbol() << "\t";
|
---|
338 | }
|
---|
339 | break;
|
---|
340 | case TremoloKey::Id :
|
---|
341 | LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getId()+1);
|
---|
342 | *file << getLocalId(currentAtom->getId()) << "\t";
|
---|
343 | break;
|
---|
344 | case TremoloKey::neighbors :
|
---|
345 | LOG(3, "Writing type " << knownKeyNames[currentField]);
|
---|
346 | writeNeighbors(file, atoi(it->substr(it->find("=") + 1, 1).c_str()), currentAtom);
|
---|
347 | break;
|
---|
348 | case TremoloKey::resSeq :
|
---|
349 | if (additionalAtomData.count(currentAtom->getId())) {
|
---|
350 | LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
|
---|
351 | *file << additionalAtomData[currentAtom->getId()].get(currentField);
|
---|
352 | } else if (currentAtom->getMolecule() != NULL) {
|
---|
353 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " its own id: " << currentAtom->getMolecule()->getId()+1);
|
---|
354 | *file << setw(4) << currentAtom->getMolecule()->getId()+1;
|
---|
355 | } else {
|
---|
356 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value: " << defaultAdditionalData.get(currentField));
|
---|
357 | *file << defaultAdditionalData.get(currentField);
|
---|
358 | }
|
---|
359 | *file << "\t";
|
---|
360 | break;
|
---|
361 | case TremoloKey::charge :
|
---|
362 | if (currentAtom->getCharge() == 0.) {
|
---|
363 | if (additionalAtomData.count(currentAtom->getId())) {
|
---|
364 | LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
|
---|
365 | *file << additionalAtomData[currentAtom->getId()].get(currentField);
|
---|
366 | } else if (additionalAtomData.count(currentAtom->GetTrueFather()->getId())) {
|
---|
367 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " stuff from father: " << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField));
|
---|
368 | *file << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField);
|
---|
369 | } else {
|
---|
370 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " AtomInfo::charge : " << currentAtom->getCharge());
|
---|
371 | *file << currentAtom->getCharge();
|
---|
372 | }
|
---|
373 | } else {
|
---|
374 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " AtomInfo::charge : " << currentAtom->getCharge());
|
---|
375 | *file << currentAtom->getCharge();
|
---|
376 | }
|
---|
377 | *file << "\t";
|
---|
378 | break;
|
---|
379 | default :
|
---|
380 | if (additionalAtomData.count(currentAtom->getId())) {
|
---|
381 | LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
|
---|
382 | *file << additionalAtomData[currentAtom->getId()].get(currentField);
|
---|
383 | } else if (additionalAtomData.count(currentAtom->GetTrueFather()->getId())) {
|
---|
384 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " stuff from father: " << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField));
|
---|
385 | *file << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField);
|
---|
386 | } else {
|
---|
387 | LOG(3, "Writing for type " << knownKeyNames[currentField] << " the default: " << defaultAdditionalData.get(currentField));
|
---|
388 | *file << defaultAdditionalData.get(currentField);
|
---|
389 | }
|
---|
390 | *file << "\t";
|
---|
391 | break;
|
---|
392 | }
|
---|
393 | }
|
---|
394 |
|
---|
395 | *file << std::endl;
|
---|
396 | }
|
---|
397 |
|
---|
398 | /**
|
---|
399 | * Writes the neighbor information of one atom to the provided stream.
|
---|
400 | *
|
---|
401 | * Note that ListOfBonds of WorldTime::CurrentTime is used.
|
---|
402 | *
|
---|
403 | * \param stream where to write neighbor information to
|
---|
404 | * \param number of neighbors
|
---|
405 | * \param reference to the atom of which to take the neighbor information
|
---|
406 | */
|
---|
407 | void FormatParser< tremolo >::writeNeighbors(std::ostream* file, const int numberOfNeighbors, const atom* currentAtom) {
|
---|
408 | const BondList& ListOfBonds = currentAtom->getListOfBonds();
|
---|
409 | // sort bonded indices
|
---|
410 | typedef std::set<atomId_t> sortedIndices;
|
---|
411 | sortedIndices sortedBonds;
|
---|
412 | for (BondList::const_iterator iter = ListOfBonds.begin();
|
---|
413 | iter != ListOfBonds.end(); ++iter)
|
---|
414 | sortedBonds.insert(getLocalId((*iter)->GetOtherAtom(currentAtom)->getId()));
|
---|
415 | // print indices
|
---|
416 | sortedIndices::const_iterator currentBond = sortedBonds.begin();
|
---|
417 | for (int i = 0; i < numberOfNeighbors; i++) {
|
---|
418 | *file << (currentBond != sortedBonds.end() ? (*currentBond) : 0) << "\t";
|
---|
419 | if (currentBond != sortedBonds.end())
|
---|
420 | ++currentBond;
|
---|
421 | }
|
---|
422 | }
|
---|
423 |
|
---|
424 | /**
|
---|
425 | * Stores keys from the ATOMDATA line in \a fields.
|
---|
426 | *
|
---|
427 | * \param line to parse the keys from
|
---|
428 | * \param offset with which offset the keys begin within the line
|
---|
429 | * \param fields which usedFields to use
|
---|
430 | */
|
---|
431 | void FormatParser< tremolo >::parseAtomDataKeysLine(
|
---|
432 | const std::string &line,
|
---|
433 | const int offset,
|
---|
434 | usedFields_t &fields) {
|
---|
435 | std::string keyword;
|
---|
436 | std::stringstream lineStream;
|
---|
437 |
|
---|
438 | lineStream << line.substr(offset);
|
---|
439 | while (lineStream.good()) {
|
---|
440 | lineStream >> keyword;
|
---|
441 | if (knownKeys[keyword.substr(0, keyword.find("="))] == TremoloKey::noKey) {
|
---|
442 | // TODO: throw exception about unknown key
|
---|
443 | cout << "Unknown key: " << keyword << " is not part of the tremolo format specification." << endl;
|
---|
444 | break;
|
---|
445 | }
|
---|
446 | fields.push_back(keyword);
|
---|
447 | }
|
---|
448 | //LOG(1, "INFO: " << fields);
|
---|
449 | }
|
---|
450 |
|
---|
451 | /** Sets the properties per atom to print to .data file by parsing line from
|
---|
452 | * \a atomdata_string.
|
---|
453 | *
|
---|
454 | * We just call \sa FormatParser< tremolo >::parseAtomDataKeysLine(), however
|
---|
455 | * we clear FormatParser< tremolo >::usedFields_save.
|
---|
456 | *
|
---|
457 | * @param atomdata_string line to parse with space-separated values
|
---|
458 | */
|
---|
459 | void FormatParser< tremolo >::setAtomData(const std::string &atomdata_string)
|
---|
460 | {
|
---|
461 | usedFields_save.clear();
|
---|
462 | parseAtomDataKeysLine(atomdata_string, 0, usedFields_save);
|
---|
463 | }
|
---|
464 |
|
---|
465 |
|
---|
466 | /**
|
---|
467 | * Reads one data line of a tremolo file and interprets it according to the keys
|
---|
468 | * obtained from the ATOMDATA line.
|
---|
469 | *
|
---|
470 | * \param line to parse as an atom
|
---|
471 | * \param *newmol molecule to add atom to
|
---|
472 | */
|
---|
473 | void FormatParser< tremolo >::readAtomDataLine(const std::string &line, molecule *newmol = NULL) {
|
---|
474 | std::stringstream lineStream;
|
---|
475 | atom* newAtom = World::getInstance().createAtom();
|
---|
476 | const atomId_t atomid = newAtom->getId();
|
---|
477 | additionalAtomData[atomid] = TremoloAtomInfoContainer(); // fill with default values
|
---|
478 | TremoloAtomInfoContainer *atomInfo = &additionalAtomData[atomid];
|
---|
479 | TremoloKey::atomDataKey currentField;
|
---|
480 | ConvertTo<double> toDouble;
|
---|
481 | ConvertTo<int> toInt;
|
---|
482 | Vector tempVector;
|
---|
483 |
|
---|
484 | // setup tokenizer, splitting up white-spaced entries
|
---|
485 | typedef boost::tokenizer<boost::char_separator<char> >
|
---|
486 | tokenizer;
|
---|
487 | boost::char_separator<char> whitespacesep(" \t");
|
---|
488 | tokenizer tokens(line, whitespacesep);
|
---|
489 | ASSERT(tokens.begin() != tokens.end(),
|
---|
490 | "FormatParser< tremolo >::readAtomDataLine - empty string, need at least ' '!");
|
---|
491 | tokenizer::iterator tok_iter = tokens.begin();
|
---|
492 | // then associate each token to each file
|
---|
493 | for (usedFields_t::const_iterator it = usedFields_load.begin(); it < usedFields_load.end(); it++) {
|
---|
494 | const std::string keyName = it->substr(0, it->find("="));
|
---|
495 | currentField = knownKeys[keyName];
|
---|
496 | const std::string word = *tok_iter;
|
---|
497 | LOG(4, "INFO: Parsing key " << keyName << " with remaining data " << word);
|
---|
498 | switch (currentField) {
|
---|
499 | case TremoloKey::x :
|
---|
500 | // for the moment, assume there are always three dimensions
|
---|
501 | for (int i=0;i<NDIM;i++) {
|
---|
502 | ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for x["+toString(i)+"]!");
|
---|
503 | LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
|
---|
504 | newAtom->set(i, toDouble(*tok_iter));
|
---|
505 | tok_iter++;
|
---|
506 | }
|
---|
507 | break;
|
---|
508 | case TremoloKey::u :
|
---|
509 | // for the moment, assume there are always three dimensions
|
---|
510 | for (int i=0;i<NDIM;i++) {
|
---|
511 | ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for u["+toString(i)+"]!");
|
---|
512 | LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
|
---|
513 | tempVector[i] = toDouble(*tok_iter);
|
---|
514 | tok_iter++;
|
---|
515 | }
|
---|
516 | newAtom->setAtomicVelocity(tempVector);
|
---|
517 | break;
|
---|
518 | case TremoloKey::type :
|
---|
519 | {
|
---|
520 | ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
|
---|
521 | LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
|
---|
522 | std::string element;
|
---|
523 | try {
|
---|
524 | element = knownTypes.getType(*tok_iter);
|
---|
525 | } catch(IllegalParserKeyException) {
|
---|
526 | // clean up
|
---|
527 | World::getInstance().destroyAtom(newAtom);
|
---|
528 | // give an error
|
---|
529 | ELOG(0, "TremoloParser: I do not understand the element token " << *tok_iter << ".");
|
---|
530 | }
|
---|
531 | // put type name into container for later use
|
---|
532 | atomInfo->set(currentField, *tok_iter);
|
---|
533 | LOG(4, "INFO: Parsing element " << (*tok_iter) << " as " << element << " according to KnownTypes.");
|
---|
534 | tok_iter++;
|
---|
535 | newAtom->setType(World::getInstance().getPeriode()->FindElement(element));
|
---|
536 | ASSERT(newAtom->getType(), "Type was not set for this atom");
|
---|
537 | break;
|
---|
538 | }
|
---|
539 | case TremoloKey::Id :
|
---|
540 | ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
|
---|
541 | LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
|
---|
542 | associateLocaltoGlobalId(toInt(*tok_iter), atomid);
|
---|
543 | tok_iter++;
|
---|
544 | break;
|
---|
545 | case TremoloKey::neighbors :
|
---|
546 | for (int i=0;i<atoi(it->substr(it->find("=") + 1, 1).c_str());i++) {
|
---|
547 | ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
|
---|
548 | LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
|
---|
549 | lineStream << *tok_iter << "\t";
|
---|
550 | tok_iter++;
|
---|
551 | }
|
---|
552 | readNeighbors(&lineStream,
|
---|
553 | atoi(it->substr(it->find("=") + 1, 1).c_str()), atomid);
|
---|
554 | break;
|
---|
555 | case TremoloKey::charge :
|
---|
556 | ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
|
---|
557 | LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
|
---|
558 | atomInfo->set(currentField, *tok_iter);
|
---|
559 | newAtom->setCharge(boost::lexical_cast<double>(*tok_iter));
|
---|
560 | tok_iter++;
|
---|
561 | break;
|
---|
562 | default :
|
---|
563 | ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
|
---|
564 | LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
|
---|
565 | atomInfo->set(currentField, *tok_iter);
|
---|
566 | tok_iter++;
|
---|
567 | break;
|
---|
568 | }
|
---|
569 | }
|
---|
570 | LOG(3, "INFO: Parsed atom " << atomid << ".");
|
---|
571 | if (newmol != NULL)
|
---|
572 | newmol->AddAtom(newAtom);
|
---|
573 | }
|
---|
574 |
|
---|
575 | bool FormatParser< tremolo >::saveAtomsInExttypes(std::ostream &output, const std::vector<atom*> &atoms, const int id) const
|
---|
576 | {
|
---|
577 | bool status = true;
|
---|
578 | // parse the file
|
---|
579 | for (std::vector<atom *>::const_iterator iter = atoms.begin();
|
---|
580 | iter != atoms.end(); ++iter) {
|
---|
581 | const int atomicid = getLocalId((*iter)->getId());
|
---|
582 | if (atomicid == -1)
|
---|
583 | status = false;
|
---|
584 | output << atomicid << "\t" << id << std::endl;
|
---|
585 | }
|
---|
586 |
|
---|
587 | return status;
|
---|
588 | }
|
---|
589 |
|
---|
590 | /**
|
---|
591 | * Reads neighbor information for one atom from the input.
|
---|
592 | *
|
---|
593 | * \param line stream where to read the information from
|
---|
594 | * \param numberOfNeighbors number of neighbors to read
|
---|
595 | * \param atomid world id of the atom the information belongs to
|
---|
596 | */
|
---|
597 | void FormatParser< tremolo >::readNeighbors(std::stringstream* line, const int numberOfNeighbors, const int atomId) {
|
---|
598 | int neighborId = 0;
|
---|
599 | for (int i = 0; i < numberOfNeighbors; i++) {
|
---|
600 | *line >> neighborId;
|
---|
601 | // 0 is used to fill empty neighbor positions in the tremolo file.
|
---|
602 | if (neighborId > 0) {
|
---|
603 | LOG(4, "INFO: Atom with global id " << atomId
|
---|
604 | << " has neighbour with serial " << neighborId);
|
---|
605 | additionalAtomData[atomId].neighbors.push_back(neighborId);
|
---|
606 | }
|
---|
607 | }
|
---|
608 | }
|
---|
609 |
|
---|
610 | /**
|
---|
611 | * Checks whether the provided name is within \a fields.
|
---|
612 | *
|
---|
613 | * \param fields which usedFields to use
|
---|
614 | * \param fieldName name to check
|
---|
615 | * \return true if the field name is used
|
---|
616 | */
|
---|
617 | bool FormatParser< tremolo >::isUsedField(const usedFields_t &fields, const std::string &fieldName) const
|
---|
618 | {
|
---|
619 | bool fieldNameExists = false;
|
---|
620 | for (usedFields_t::const_iterator usedField = fields.begin();
|
---|
621 | usedField != fields.end(); usedField++) {
|
---|
622 | if (usedField->substr(0, usedField->find("=")) == fieldName)
|
---|
623 | fieldNameExists = true;
|
---|
624 | }
|
---|
625 |
|
---|
626 | return fieldNameExists;
|
---|
627 | }
|
---|
628 |
|
---|
629 |
|
---|
630 | /**
|
---|
631 | * Adds the collected neighbor information to the atoms in the world. The atoms
|
---|
632 | * are found by their current ID and mapped to the corresponding atoms with the
|
---|
633 | * Id found in the parsed file.
|
---|
634 | *
|
---|
635 | * @param atoms vector with all newly added (global) atomic ids
|
---|
636 | */
|
---|
637 | void FormatParser< tremolo >::processNeighborInformation(const std::vector<atomId_t> &atoms) {
|
---|
638 | if (!isUsedField(usedFields_load, "neighbors")) {
|
---|
639 | return;
|
---|
640 | }
|
---|
641 |
|
---|
642 | for (std::vector<atomId_t>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
|
---|
643 | ASSERT(additionalAtomData.count(*iter) != 0,
|
---|
644 | "FormatParser< tremolo >::processNeighborInformation() - global id "
|
---|
645 | +toString(*iter)+" unknown in additionalAtomData.");
|
---|
646 | TremoloAtomInfoContainer ¤tInfo = additionalAtomData[*iter];
|
---|
647 | ASSERT (!currentInfo.neighbors_processed,
|
---|
648 | "FormatParser< tremolo >::processNeighborInformation() - neighbors of new atom "
|
---|
649 | +toString(*iter)+" are already processed.");
|
---|
650 | for(std::vector<int>::const_iterator neighbor = currentInfo.neighbors.begin();
|
---|
651 | neighbor != currentInfo.neighbors.end(); neighbor++
|
---|
652 | ) {
|
---|
653 | LOG(3, "INFO: Creating bond between ("
|
---|
654 | << *iter
|
---|
655 | << ") and ("
|
---|
656 | << getGlobalId(*neighbor) << "|" << *neighbor << ")");
|
---|
657 | ASSERT(getGlobalId(*neighbor) != -1,
|
---|
658 | "FormatParser< tremolo >::processNeighborInformation() - global id to local id "
|
---|
659 | +toString(*neighbor)+" is unknown.");
|
---|
660 | World::getInstance().getAtom(AtomById(*iter))
|
---|
661 | ->addBond(WorldTime::getTime(), World::getInstance().getAtom(AtomById(getGlobalId(*neighbor))));
|
---|
662 | }
|
---|
663 | currentInfo.neighbors_processed = true;
|
---|
664 | }
|
---|
665 | }
|
---|
666 |
|
---|
667 | /**
|
---|
668 | * Replaces atom IDs read from the file by the corresponding world IDs. All IDs
|
---|
669 | * IDs of the input string will be replaced; expected separating characters are
|
---|
670 | * "-" and ",".
|
---|
671 | *
|
---|
672 | * \param string in which atom IDs should be adapted
|
---|
673 | *
|
---|
674 | * \return input string with modified atom IDs
|
---|
675 | */
|
---|
676 | std::string FormatParser< tremolo >::adaptIdDependentDataString(std::string data) {
|
---|
677 | // there might be no IDs
|
---|
678 | if (data == "-") {
|
---|
679 | return "-";
|
---|
680 | }
|
---|
681 |
|
---|
682 | char separator;
|
---|
683 | int id;
|
---|
684 | std::stringstream line, result;
|
---|
685 | line << data;
|
---|
686 |
|
---|
687 | line >> id;
|
---|
688 | result << getGlobalId(id);
|
---|
689 | while (line.good()) {
|
---|
690 | line >> separator >> id;
|
---|
691 | result << separator << getGlobalId(id);
|
---|
692 | }
|
---|
693 |
|
---|
694 | return result.str();
|
---|
695 | }
|
---|
696 |
|
---|
697 | /**
|
---|
698 | * Corrects the atom IDs in each imprData entry to the corresponding world IDs
|
---|
699 | * as they might differ from the originally read IDs.
|
---|
700 | */
|
---|
701 | void FormatParser< tremolo >::adaptImprData() {
|
---|
702 | if (!isUsedField(usedFields_load, "imprData")) {
|
---|
703 | return;
|
---|
704 | }
|
---|
705 |
|
---|
706 | for(std::map<const atomId_t, TremoloAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
|
---|
707 | currentInfo != additionalAtomData.end(); currentInfo++
|
---|
708 | ) {
|
---|
709 | currentInfo->second.imprData = adaptIdDependentDataString(currentInfo->second.imprData);
|
---|
710 | }
|
---|
711 | }
|
---|
712 |
|
---|
713 | /**
|
---|
714 | * Corrects the atom IDs in each torsion entry to the corresponding world IDs
|
---|
715 | * as they might differ from the originally read IDs.
|
---|
716 | */
|
---|
717 | void FormatParser< tremolo >::adaptTorsion() {
|
---|
718 | if (!isUsedField(usedFields_load, "torsion")) {
|
---|
719 | return;
|
---|
720 | }
|
---|
721 |
|
---|
722 | for(std::map<const atomId_t, TremoloAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
|
---|
723 | currentInfo != additionalAtomData.end(); currentInfo++
|
---|
724 | ) {
|
---|
725 | currentInfo->second.torsion = adaptIdDependentDataString(currentInfo->second.torsion);
|
---|
726 | }
|
---|
727 | }
|
---|
728 |
|
---|