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