source: src/Parser/TremoloParser.cpp@ cd91bd

Candidate_v1.6.1 ChemicalSpaceEvaluator PythonUI_with_named_parameters TremoloParser_IncreasedPrecision
Last change on this file since cd91bd was 33af20, checked in by Frederik Heber <frederik.heber@…>, 8 years ago

TremoloParser::save() now supports multiple time steps.

  • TESTS: Removed XFAIL from tests Tremolo with multiple time steps.
  • Property mode set to 100644
File size: 39.5 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/*
25 * TremoloParser.cpp
26 *
27 * Created on: Mar 2, 2010
28 * Author: metzler
29 */
30
31// include config.h
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36//#include "CodePatterns/MemDebug.hpp"
37
38#include "CodePatterns/Assert.hpp"
39#include "CodePatterns/Log.hpp"
40#include "CodePatterns/toString.hpp"
41#include "CodePatterns/Verbose.hpp"
42
43#include "TremoloParser.hpp"
44
45#include "Atom/atom.hpp"
46#include "Bond/bond.hpp"
47#include "Box.hpp"
48#include "Descriptors/AtomIdDescriptor.hpp"
49#include "Element/element.hpp"
50#include "Element/periodentafel.hpp"
51#include "Parser/Exceptions.hpp"
52#include "Potentials/Particles/ParticleRegistry.hpp"
53#include "LinearAlgebra/RealSpaceMatrix.hpp"
54#include "Potentials/Particles/ParticleRegistry.hpp"
55#include "molecule.hpp"
56#include "World.hpp"
57#include "WorldTime.hpp"
58
59
60#include <algorithm>
61#include <boost/bind.hpp>
62#include <boost/function.hpp>
63#include <boost/lambda/lambda.hpp>
64#include <boost/lexical_cast.hpp>
65#include <boost/tokenizer.hpp>
66#include <iostream>
67#include <iomanip>
68#include <map>
69#include <sstream>
70#include <string>
71#include <vector>
72
73#include <boost/assign/list_of.hpp> // for 'map_list_of()'
74#include <boost/assert.hpp>
75
76// declare specialized static variables
77const std::string FormatParserTrait<tremolo>::name = "tremolo";
78const std::string FormatParserTrait<tremolo>::suffix = "data";
79const ParserTypes FormatParserTrait<tremolo>::type = tremolo;
80
81// static instances
82std::map<std::string, TremoloKey::atomDataKey> FormatParser<tremolo>::knownKeys =
83 boost::assign::map_list_of("x",TremoloKey::x)
84 ("u",TremoloKey::u)
85 ("F",TremoloKey::F)
86 ("stress",TremoloKey::stress)
87 ("Id",TremoloKey::Id)
88 ("neighbors",TremoloKey::neighbors)
89 ("imprData",TremoloKey::imprData)
90 ("GroupMeasureTypeNo",TremoloKey::GroupMeasureTypeNo)
91 ("type",TremoloKey::type)
92 ("extType",TremoloKey::extType)
93 ("name",TremoloKey::name)
94 ("resName",TremoloKey::resName)
95 ("chainID",TremoloKey::chainID)
96 ("resSeq",TremoloKey::resSeq)
97 ("occupancy",TremoloKey::occupancy)
98 ("tempFactor",TremoloKey::tempFactor)
99 ("segID",TremoloKey::segID)
100 ("Charge",TremoloKey::Charge)
101 ("charge",TremoloKey::charge)
102 ("GrpTypeNo",TremoloKey::GrpTypeNo)
103 ("torsion",TremoloKey::torsion)
104 (" ",TremoloKey::noKey); // with this we can detect invalid keys
105
106/**
107 * Constructor.
108 */
109FormatParser< tremolo >::FormatParser() :
110 FormatParser_common(NULL),
111 idglobalizer(boost::bind(&FormatParser< tremolo >::getGlobalId, this, _1)),
112 idlocalizer(boost::bind(&FormatParser< tremolo >::getLocalId, this, _1))
113{
114 // invert knownKeys for debug output
115 for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin(); iter != knownKeys.end(); ++iter)
116 knownKeyNames.insert( make_pair( iter->second, iter->first) );
117
118 additionalAtomData.clear();
119}
120
121
122/**
123 * Destructor.
124 */
125FormatParser< tremolo >::~FormatParser()
126{
127 usedFields_save.clear();
128 additionalAtomData.clear();
129}
130
131/**
132 * Loads atoms from a tremolo-formatted file.
133 *
134 * \param tremolo file
135 */
136void FormatParser< tremolo >::load(istream* file) {
137 std::string line;
138 std::string::size_type location;
139
140 // reset the id maps
141 resetIdAssociations();
142
143 molecule *newmol = World::getInstance().createMolecule();
144 newmol->ActiveFlag = true;
145 usedFields_t usedFields_temp;
146 size_t timestep = 0;
147 atom *addedatom = NULL;
148 std::list<atom *> addedatoms;
149 std::list<atom *>::iterator atomiter = addedatoms.begin();
150 while (file->good()) {
151 std::getline(*file, line, '\n');
152 // we only parse in the first ATOMDATA line
153 location = line.find("ATOMDATA", 0);
154 if (location != string::npos) {
155 parseAtomDataKeysLine(line, location + 8, usedFields_temp);
156 if (usedFields_load.empty()) {
157 // first ATOMDATA: use this line
158 usedFields_load.insert(usedFields_load.end(), usedFields_temp.begin(), usedFields_temp.end());
159 LOG(3, "DEBUG: Local usedFields is: " << usedFields_load);
160 } else {
161 // following ATOMDATA: check against present
162 LOG(3, "DEBUG: Parsed check usedFields is: " << usedFields_temp);
163 if (usedFields_load != usedFields_temp) {
164 ELOG(1, "File contains multiple time steps with differing ATOMDATA line, parsing only first time step.");
165 break;
166 } else {
167 // update time step
168 ++timestep;
169 atomiter = addedatoms.begin();
170 }
171 }
172 usedFields_temp.clear();
173 } else {
174 if (line.length() > 0 && line.at(0) != '#') {
175 if (timestep != 0) {
176 ASSERT( atomiter != addedatoms.end(),
177 "FormatParser< tremolo >::load() - timesteps contains more atoms than first one.");
178 addedatom = *atomiter++;
179 }
180 readAtomDataLine(line, newmol, timestep, addedatom);
181 if (timestep == 0)
182 addedatoms.push_back(addedatom);
183 }
184 }
185 }
186 ASSERT( addedatoms.size() == newmol->size(),
187 "FormatParser< tremolo >::load() - number of atoms in mol and addedatoms differ.");
188
189 // refresh atom::nr and atom::name
190 std::vector<atomId_t> atoms(newmol->getAtomCount());
191 std::transform(newmol->begin(), newmol->end(), atoms.begin(),
192 boost::bind(&atom::getId, _1));
193 processNeighborInformation(atoms);
194 adaptImprData(atoms);
195 adaptTorsion(atoms);
196
197 // append usedFields to global usedFields, is made unique on save, clear after use
198 usedFields_save.insert(usedFields_save.end(), usedFields_load.begin(), usedFields_load.end());
199 usedFields_load.clear();
200}
201
202/**
203 * Saves the \a atoms into as a tremolo file.
204 *
205 * \param file where to save the state
206 * \param atoms atoms to store
207 */
208void FormatParser< tremolo >::save(
209 std::ostream* file,
210 const std::vector<const atom *> &AtomList) {
211 LOG(2, "DEBUG: Saving changes to tremolo.");
212
213 // install default usedFields if empty so far
214 if (usedFields_save.empty()) {
215 // default behavior: use all possible keys on output
216 for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin();
217 iter != knownKeys.end(); ++iter)
218 if (iter->second != TremoloKey::noKey) // don't add noKey
219 usedFields_save.push_back(iter->first);
220 }
221 // make present usedFields_save unique
222 makeUsedFieldsUnique(usedFields_save);
223 LOG(1, "INFO: Global (with unique entries) usedFields_save is: " << usedFields_save);
224
225 // distribute ids continuously
226 distributeContinuousIds(AtomList);
227
228 std::pair<size_t, size_t> minmax_trajectories =
229 getMinMaxTrajectories(AtomList);
230 LOG(2, "INFO: There are " << minmax_trajectories.second << " steps to save.");
231
232 for (size_t step = 0; (step < minmax_trajectories.second) || (step == 0); ++step) {
233 // store atomdata
234 save_AtomDataLine(file);
235
236 // store box only on first step
237 if (step == 0)
238 save_BoxLine(file);
239
240 // store particles
241 for (std::vector<const atom*>::const_iterator atomIt = AtomList.begin();
242 atomIt != AtomList.end(); ++atomIt)
243 saveLine(file, *atomIt, step);
244 }
245}
246
247struct usedFieldsWeakComparator
248{
249 /** Special comparator regards "neighbors=4" and "neighbors=2" as equal
250 *
251 * \note This one is used for making usedFields unique, i.e. throwing out the "smaller"
252 * neighbors.
253 */
254 bool operator()(const std::string &a, const std::string &b) const
255 {
256 // only compare up to first equality sign
257 return (a.substr(0, a.find_first_of('=')) == b.substr(0, b.find_first_of('=')));
258 }
259};
260
261struct usedFieldsSpecialOrderer
262{
263 /** Special string comparator that regards "neighbors=4" < "neighbors=2" as true and
264 * the other way round as false.
265 *
266 * Here, we implement the operator "\a < \b" in a special way to allow the
267 * above.
268 *
269 * \note This one is used for sorting usedFields in preparation for making it unique.
270 */
271 bool operator()(const std::string &a, const std::string &b) const
272 {
273 // only compare up to first equality sign
274 size_t a_equality = a.find_first_of('=');
275 size_t b_equality = b.find_first_of('=');
276 // if key before equality is not equal, return whether it is smaller or not
277 if (a.substr(0, a_equality) != b.substr(0, b_equality)) {
278 return a.substr(0, a_equality) < b.substr(0, b_equality);
279 } else { // now we know that the key before equality is the same in either string
280 // if one of them has no equality, the one with equality must go before
281 if ((a_equality != std::string::npos) && (b_equality == std::string::npos))
282 return true;
283 if ((a_equality == std::string::npos) && (b_equality != std::string::npos))
284 return false;
285 // if both don't have equality (and the token before is equal), it is not "<" but "=="
286 if ((a_equality == std::string::npos) && (b_equality == std::string::npos))
287 return false;
288 // if now both have equality sign, the larger value after it, must come first
289 return a.substr(a_equality, std::string::npos) > b.substr(b_equality, std::string::npos);
290 }
291 }
292};
293
294/** Helper function to make \given fields unique while preserving the order of first appearance.
295 *
296 * As std::unique only removes element if equal to predecessor, a vector is only
297 * made unique if sorted beforehand. But sorting would destroy order of first
298 * appearance, hence we do the sorting on a temporary field and add the unique
299 * elements in the order as in \a fields.
300 *
301 * @param fields usedFields to make unique while preserving order of appearance
302 */
303void FormatParser< tremolo >::makeUsedFieldsUnique(usedFields_t &fields) const
304{
305 // std::unique only removes if predecessor is equal, not over whole range, hence do it manually
306 usedFields_t temp_fields(fields);
307 usedFieldsSpecialOrderer SpecialOrderer;
308 usedFieldsWeakComparator WeakComparator;
309 std::sort(temp_fields.begin(), temp_fields.end(), SpecialOrderer);
310 usedFields_t::iterator it =
311 std::unique(temp_fields.begin(), temp_fields.end(), WeakComparator);
312 temp_fields.erase(it, temp_fields.end());
313 usedFields_t usedfields(fields);
314 fields.clear();
315 fields.reserve(temp_fields.size());
316 // now go through each usedFields entry, check if in temp_fields and remove there on first occurence
317 for (usedFields_t::const_iterator iter = usedfields.begin();
318 iter != usedfields.end(); ++iter) {
319 usedFields_t::iterator uniqueiter =
320 std::find(temp_fields.begin(), temp_fields.end(), *iter);
321 if (uniqueiter != temp_fields.end()) {
322 fields.push_back(*iter);
323 // add only once to ATOMDATA
324 temp_fields.erase(uniqueiter);
325 }
326 }
327 ASSERT( temp_fields.empty(),
328 "FormatParser< tremolo >::save() - still unique entries left in temp_fields after unique?");
329}
330
331
332/** Resets and distributes the indices continuously.
333 *
334 * \param atoms atoms to store
335 */
336void FormatParser< tremolo >::distributeContinuousIds(
337 const std::vector<const atom *> &AtomList)
338{
339 resetIdAssociations();
340 atomId_t lastid = 0;
341 for (std::vector<const atom*>::const_iterator atomIt = AtomList.begin();
342 atomIt != AtomList.end(); ++atomIt)
343 associateLocaltoGlobalId(++lastid, (*atomIt)->getId());
344}
345
346/** Store Atomdata line to \a file.
347 *
348 * @param file output stream
349 */
350void FormatParser< tremolo >::save_AtomDataLine(std::ostream* file) const
351{
352 *file << "# ATOMDATA";
353 for (usedFields_t::const_iterator it=usedFields_save.begin();
354 it != usedFields_save.end(); ++it)
355 *file << "\t" << *it;
356 *file << std::endl;
357}
358
359/** Store Box info to \a file
360 *
361 * @param file output stream
362 */
363void FormatParser< tremolo >::save_BoxLine(std::ostream* file) const
364{
365 *file << "# Box";
366 const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
367 for (size_t i=0; i<NDIM;++i)
368 for (size_t j=0; j<NDIM;++j)
369 *file << "\t" << M.at(i,j);
370 *file << std::endl;
371}
372
373/** Add default info, when new atom is added to World.
374 *
375 * @param id of atom
376 */
377void FormatParser< tremolo >::AtomInserted(atomId_t id)
378{
379 std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
380 ASSERT(iter == additionalAtomData.end(),
381 "FormatParser< tremolo >::AtomInserted() - additionalAtomData already present for newly added atom "
382 +toString(id)+".");
383 // don't add entry, as this gives a default resSeq of 0 not the molecule id
384 // additionalAtomData.insert( std::make_pair(id, TremoloAtomInfoContainer()) );
385}
386
387/** Remove additional AtomData info, when atom has been removed from World.
388 *
389 * @param id of atom
390 */
391void FormatParser< tremolo >::AtomRemoved(atomId_t id)
392{
393 std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
394 // as we do not insert AtomData on AtomInserted, we cannot be assured of its presence
395// ASSERT(iter != additionalAtomData.end(),
396// "FormatParser< tremolo >::AtomRemoved() - additionalAtomData is not present for atom "
397// +toString(id)+" to remove.");
398 if (iter != additionalAtomData.end())
399 additionalAtomData.erase(iter);
400}
401
402template <class T>
403T NoOp(const atom * const)
404{
405 return T();
406}
407
408template <class T>
409void writeEntryFromAdditionalAtomData_ifpresent(
410 std::map<const atomId_t, TremoloAtomInfoContainer> &_additionalAtomData,
411 std::map<TremoloKey::atomDataKey, std::string> &_knownKeyNames,
412 std::ostream* _file,
413 const atom * const _currentAtom,
414 const TremoloKey::atomDataKey _currentField,
415 const typename boost::function<T (const atom * const)> _getter)
416{
417 if (_additionalAtomData.count(_currentAtom->getId())) {
418 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << ": "
419 << _additionalAtomData[_currentAtom->getId()].get(_currentField));
420 *_file << _additionalAtomData[_currentAtom->getId()].get(_currentField);
421 } else if (_additionalAtomData.count(_currentAtom->GetTrueFather()->getId())) {
422 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " stuff from father: "
423 << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField));
424 *_file << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField);
425 } else {
426 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " its default value: "
427 << _getter(_currentAtom));
428 *_file << _getter(_currentAtom);
429 }
430 *_file << "\t";
431}
432
433template <class T>
434void writeAdditionalAtomDataEntry(
435 std::map<const atomId_t, TremoloAtomInfoContainer> &_additionalAtomData,
436 std::map<TremoloKey::atomDataKey, std::string> &_knownKeyNames,
437 TremoloAtomInfoContainer &_defaultAdditionalData,
438 std::ostream* _file,
439 const atom * const _currentAtom,
440 const TremoloKey::atomDataKey _currentField)
441{
442 if (_additionalAtomData.count(_currentAtom->getId())) {
443 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << ": "
444 << _additionalAtomData[_currentAtom->getId()].get(_currentField));
445 *_file << _additionalAtomData[_currentAtom->getId()].get(_currentField);
446 } else if (_additionalAtomData.count(_currentAtom->GetTrueFather()->getId())) {
447 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " stuff from father: "
448 << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField));
449 *_file << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField);
450 } else {
451 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " its default value: "
452 << _defaultAdditionalData.get(_currentField));
453 *_file << _defaultAdditionalData.get(_currentField);
454 }
455 *_file << "\t";
456}
457
458template <class T>
459void writeEntryFromAdditionalAtomData_ifnotempty(
460 std::map<const atomId_t, TremoloAtomInfoContainer> &_additionalAtomData,
461 std::map<TremoloKey::atomDataKey, std::string> &_knownKeyNames,
462 TremoloAtomInfoContainer &_defaultAdditionalData,
463 std::ostream* _file,
464 const atom * const _currentAtom,
465 const TremoloKey::atomDataKey _currentField,
466 const typename boost::function<T (const atom * const)> _getter)
467{
468 if (_additionalAtomData.count(_currentAtom->getId())) {
469 if (_additionalAtomData[_currentAtom->getId()].get(_currentField) != "-") {
470 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << ": "
471 << _additionalAtomData[_currentAtom->getId()].get(_currentField));
472 *_file << _additionalAtomData[_currentAtom->getId()].get(_currentField);
473 } else {
474 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " default value: "
475 << _getter(_currentAtom));
476 *_file << _getter(_currentAtom);
477 }
478 } else if (_additionalAtomData.count(_currentAtom->GetTrueFather()->getId())) {
479 if (_additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField) != "-") {
480 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " stuff from father: "
481 << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField));
482 *_file << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField);
483 } else {
484 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " default value from father: "
485 << _getter(_currentAtom->GetTrueFather()));
486 *_file << _getter(_currentAtom->GetTrueFather());
487 }
488 } else {
489 LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " its default value: "
490 << _getter(_currentAtom));
491 *_file << _getter(_currentAtom);
492 }
493 *_file << "\t";
494}
495
496/**
497 * Writes one line of tremolo-formatted data to the provided stream.
498 *
499 * \param stream where to write the line to
500 * \param reference to the atom of which information should be written
501 */
502void FormatParser< tremolo >::saveLine(
503 std::ostream* file,
504 const atom * const currentAtom,
505 const size_t _timestep)
506{
507 TremoloKey::atomDataKey currentField;
508
509 LOG(4, "INFO: Saving atom " << *currentAtom << ", its father id is "
510 << currentAtom->GetTrueFather()->getId() << " at time step " << _timestep);
511
512 for (usedFields_t::iterator it = usedFields_save.begin(); it != usedFields_save.end(); it++) {
513 currentField = knownKeys[it->substr(0, it->find("="))];
514 switch (currentField) {
515 case TremoloKey::x :
516 // for the moment, assume there are always three dimensions
517 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getPositionAtStep(_timestep));
518 *file << currentAtom->atStep(0, _timestep) << "\t";
519 *file << currentAtom->atStep(1, _timestep) << "\t";
520 *file << currentAtom->atStep(2, _timestep) << "\t";
521 break;
522 case TremoloKey::u :
523 // for the moment, assume there are always three dimensions
524 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicVelocityAtStep(_timestep));
525 *file << currentAtom->getAtomicVelocityAtStep(_timestep)[0] << "\t";
526 *file << currentAtom->getAtomicVelocityAtStep(_timestep)[1] << "\t";
527 *file << currentAtom->getAtomicVelocityAtStep(_timestep)[2] << "\t";
528 break;
529 case TremoloKey::F :
530 // for the moment, assume there are always three dimensions
531 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicForceAtStep(_timestep));
532 *file << currentAtom->getAtomicForceAtStep(_timestep)[0] << "\t";
533 *file << currentAtom->getAtomicForceAtStep(_timestep)[1] << "\t";
534 *file << currentAtom->getAtomicForceAtStep(_timestep)[2] << "\t";
535 break;
536 case TremoloKey::type :
537 writeEntryFromAdditionalAtomData_ifnotempty<std::string>(
538 additionalAtomData,
539 knownKeyNames,
540 defaultAdditionalData,
541 file,
542 currentAtom,
543 currentField,
544 boost::bind(&element::getSymbol, boost::bind(&AtomInfo::getType, _1)));
545 break;
546 case TremoloKey::Id :
547 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getId()+1);
548 *file << getLocalId(currentAtom->getId()) << "\t";
549 break;
550 case TremoloKey::neighbors :
551 LOG(3, "Writing type " << knownKeyNames[currentField]);
552 writeNeighbors(file, atoi(it->substr(it->find("=") + 1, 1).c_str()), currentAtom);
553 break;
554 case TremoloKey::imprData :
555 case TremoloKey::torsion :
556 LOG(3, "Writing type " << knownKeyNames[currentField]);
557 *file << adaptIdDependentDataString(
558 additionalAtomData[currentAtom->getId()].get(currentField),
559 idlocalizer)
560 << "\t";
561 break;
562 case TremoloKey::resSeq :
563 if (additionalAtomData.count(currentAtom->getId())) {
564 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
565 *file << additionalAtomData[currentAtom->getId()].get(currentField);
566 } else if (currentAtom->getMolecule() != NULL) {
567 LOG(3, "Writing for type " << knownKeyNames[currentField] << " its own id: " << currentAtom->getMolecule()->getId()+1);
568 *file << setw(4) << currentAtom->getMolecule()->getId();
569 } else {
570 LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value: " << defaultAdditionalData.get(currentField));
571 *file << defaultAdditionalData.get(currentField);
572 }
573 *file << "\t";
574 break;
575 case TremoloKey::charge :
576 if (currentAtom->getCharge() == 0.) {
577 writeEntryFromAdditionalAtomData_ifpresent<double>(
578 additionalAtomData,
579 knownKeyNames,
580 file,
581 currentAtom,
582 currentField,
583 boost::bind(&AtomInfo::getCharge, _1));
584 } else {
585 LOG(3, "Writing for type " << knownKeyNames[currentField] << " AtomInfo::charge : " << currentAtom->getCharge());
586 *file << currentAtom->getCharge() << "\t";
587 }
588 break;
589 default :
590 writeAdditionalAtomDataEntry<std::string>(
591 additionalAtomData,
592 knownKeyNames,
593 defaultAdditionalData,
594 file,
595 currentAtom,
596 currentField);
597 break;
598 }
599 }
600
601 *file << std::endl;
602}
603
604/**
605 * Writes the neighbor information of one atom to the provided stream.
606 *
607 * Note that ListOfBonds of WorldTime::CurrentTime is used.
608 *
609 * \param stream where to write neighbor information to
610 * \param number of neighbors
611 * \param reference to the atom of which to take the neighbor information
612 */
613void FormatParser< tremolo >::writeNeighbors(
614 std::ostream* file,
615 const int numberOfNeighbors,
616 const atom * const currentAtom) {
617 const BondList& ListOfBonds = currentAtom->getListOfBonds();
618 // sort bonded indices
619 typedef std::set<atomId_t> sortedIndices;
620 sortedIndices sortedBonds;
621 for (BondList::const_iterator iter = ListOfBonds.begin();
622 iter != ListOfBonds.end(); ++iter)
623 sortedBonds.insert(getLocalId((*iter)->GetOtherAtom(currentAtom)->getId()));
624 // print indices
625 sortedIndices::const_iterator currentBond = sortedBonds.begin();
626 for (int i = 0; i < numberOfNeighbors; i++) {
627 *file << (currentBond != sortedBonds.end() ? (*currentBond) : 0) << "\t";
628 if (currentBond != sortedBonds.end())
629 ++currentBond;
630 }
631}
632
633/**
634 * Stores keys from the ATOMDATA line in \a fields.
635 *
636 * \param line to parse the keys from
637 * \param offset with which offset the keys begin within the line
638 * \param fields which usedFields to use
639 */
640void FormatParser< tremolo >::parseAtomDataKeysLine(
641 const std::string &line,
642 const int offset,
643 usedFields_t &fields) {
644 std::string keyword;
645 std::stringstream lineStream;
646
647 lineStream << line.substr(offset);
648 lineStream >> ws;
649 while (lineStream.good()) {
650 lineStream >> keyword;
651 if (knownKeys[keyword.substr(0, keyword.find("="))] == TremoloKey::noKey) {
652 // TODO: throw exception about unknown key
653 cout << "Unknown key: " << keyword.substr(0, keyword.find("=")) << " is not part of the tremolo format specification." << endl;
654 throw IllegalParserKeyException();
655 break;
656 }
657 fields.push_back(keyword);
658 lineStream >> ws;
659 }
660 LOG(2, "INFO: " << fields);
661}
662
663/**
664 * Tests whether the keys from the ATOMDATA line can be read correctly.
665 *
666 * \param line to parse the keys from
667 */
668bool FormatParser< tremolo >::testParseAtomDataKeysLine(
669 const std::string &line) {
670 std::string keyword;
671 std::stringstream lineStream;
672
673 // check string after ATOMDATA
674 const std::string AtomData("ATOMDATA");
675 const size_t AtomDataOffset = line.find(AtomData, 0);
676 if (AtomDataOffset == std::string::npos)
677 lineStream << line;
678 else
679 lineStream << line.substr(AtomDataOffset + AtomData.length());
680 while (lineStream.good()) {
681 lineStream >> keyword;
682 //LOG(2, "DEBUG: Checking key " << keyword.substr(0, keyword.find("=")) << ".");
683 if (knownKeys[keyword.substr(0, keyword.find("="))] == TremoloKey::noKey)
684 return false;
685 }
686 //LOG(1, "INFO: " << fields);
687 return true;
688}
689
690std::string FormatParser< tremolo >::getAtomData() const
691{
692 std::stringstream output;
693 std::for_each(usedFields_save.begin(), usedFields_save.end(),
694 output << boost::lambda::_1 << " ");
695 const std::string returnstring(output.str());
696 return returnstring.substr(0, returnstring.find_last_of(" "));
697}
698
699/** Appends the properties per atom to print to .data file by parsing line from
700 * \a atomdata_string.
701 *
702 * We just call \sa FormatParser< tremolo >::parseAtomDataKeysLine().
703 *
704 * @param atomdata_string line to parse with space-separated values
705 */
706void FormatParser< tremolo >::setAtomData(const std::string &atomdata_string)
707{
708 parseAtomDataKeysLine(atomdata_string, 0, usedFields_save);
709}
710
711/** Sets the properties per atom to print to .data file by parsing line from
712 * \a atomdata_string.
713 *
714 * We just call \sa FormatParser< tremolo >::parseAtomDataKeysLine(), however
715 * we clear FormatParser< tremolo >::usedFields_save.
716 *
717 * @param atomdata_string line to parse with space-separated values
718 */
719void FormatParser< tremolo >::resetAtomData(const std::string &atomdata_string)
720{
721 usedFields_save.clear();
722 parseAtomDataKeysLine(atomdata_string, 0, usedFields_save);
723}
724
725
726/**
727 * Reads one data line of a tremolo file and interprets it according to the keys
728 * obtained from the ATOMDATA line.
729 *
730 * \param line to parse as an atom
731 * \param *newmol molecule to add atom to
732 * \param _timestep time step to parse to
733 * \param _addedatom ref to added atom
734 */
735void FormatParser< tremolo >::readAtomDataLine(
736 const std::string &line,
737 molecule *newmol,
738 const size_t _timestep,
739 atom *&_addedatoms)
740{
741 std::stringstream lineStream;
742 atom* newAtom = NULL;
743 atomId_t atomid = -1;
744 if (_timestep == 0) {
745 newAtom = World::getInstance().createAtom();
746 _addedatoms = newAtom;
747 atomid = newAtom->getId();
748 additionalAtomData[atomid] = TremoloAtomInfoContainer(); // fill with default values
749 } else {
750 newAtom = _addedatoms;
751 atomid = newAtom->getId();
752 }
753 TremoloAtomInfoContainer *atomInfo = &additionalAtomData[atomid];
754 TremoloKey::atomDataKey currentField;
755 ConvertTo<double> toDouble;
756 ConvertTo<int> toInt;
757 Vector tempVector;
758 const ParticleRegistry& registry = ParticleRegistry::getConstInstance();
759 const periodentafel& periode = *World::getInstance().getPeriode();
760
761 // setup tokenizer, splitting up white-spaced entries
762 typedef boost::tokenizer<boost::char_separator<char> >
763 tokenizer;
764 boost::char_separator<char> whitespacesep(" \t");
765 tokenizer tokens(line, whitespacesep);
766 ASSERT(tokens.begin() != tokens.end(),
767 "FormatParser< tremolo >::readAtomDataLine - empty string, need at least ' '!");
768 tokenizer::const_iterator tok_iter = tokens.begin();
769 // then associate each token to each file
770 for (usedFields_t::const_iterator it = usedFields_load.begin(); it != usedFields_load.end(); it++) {
771 const std::string keyName = it->substr(0, it->find("="));
772 currentField = knownKeys[keyName];
773 ASSERT(tok_iter != tokens.end(),
774 "FormatParser< tremolo >::readAtomDataLine - too few entries in line '"+line+"'!");
775 const std::string &word = *tok_iter;
776 LOG(4, "INFO: Parsing key " << keyName << " with remaining data " << word);
777 switch (currentField) {
778 case TremoloKey::x :
779 // for the moment, assume there are always three dimensions
780 for (int i=0;i<NDIM;i++) {
781 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for x["+toString(i)+"]!");
782 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
783 newAtom->setAtStep(i, _timestep, toDouble(word));
784 tok_iter++;
785 }
786 break;
787 case TremoloKey::u :
788 // for the moment, assume there are always three dimensions
789 for (int i=0;i<NDIM;i++) {
790 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for u["+toString(i)+"]!");
791 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
792 tempVector[i] = toDouble(word);
793 tok_iter++;
794 }
795 newAtom->setAtomicVelocityAtStep(_timestep, tempVector);
796 break;
797 case TremoloKey::F :
798 // for the moment, assume there are always three dimensions
799 for (int i=0;i<NDIM;i++) {
800 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for F["+toString(i)+"]!");
801 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
802 tempVector[i] = toDouble(word);
803 tok_iter++;
804 }
805 newAtom->setAtomicForceAtStep(_timestep, tempVector);
806 break;
807 case TremoloKey::type :
808 {
809 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
810 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
811 std::string elementname;
812 const element * elem = NULL;
813 if (!registry.isPresentByName(word)) {
814 std::string lowercase_word(word);
815 std::transform(word.begin()+1, word.end(), lowercase_word.begin()+1, ::tolower);
816 elem = periode.FindElement(lowercase_word);
817 if (elem == NULL) {
818 // clean up
819 World::getInstance().destroyAtom(newAtom);
820 // give an error
821 ELOG(0, "TremoloParser: tokens " << word << "/" << lowercase_word
822 << " is unknown to neither ParticleRegistry nor Periodentafel.");
823 return;
824 } else {
825 elementname = elem->getSymbol();
826 }
827 } else {
828 const Particle * const p = registry.getByName(word);
829 elementname = p->getElement();
830 elem = periode.FindElement(elementname);
831 }
832 // put type name into container for later use
833 atomInfo->set(currentField, word);
834 LOG(4, "INFO: Parsing element " << (word) << " as " << elementname << " according to KnownTypes.");
835 tok_iter++;
836 newAtom->setType(elem);
837 ASSERT(newAtom->getType(), "Type was not set for this atom");
838 break;
839 }
840 case TremoloKey::Id :
841 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
842 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
843 if (_timestep == 0) {
844 associateLocaltoGlobalId(toInt(word), atomid);
845 } else {
846 // check association is the same
847 ASSERT( (atomId_t)getGlobalId(toInt(word)) == atomid,
848 "FormatParser< tremolo >::readAtomDataLine() - differing global id "+toString(atomid)
849 +" on timestep "+toString(_timestep));
850 ASSERT( getLocalId(atomid) == toInt(word),
851 "FormatParser< tremolo >::readAtomDataLine() - differing local id "+toString(toInt(word))
852 +" on timestep "+toString(_timestep));
853 }
854 tok_iter++;
855 break;
856 case TremoloKey::neighbors :
857 for (int i=0;i<atoi(it->substr(it->find("=") + 1, 1).c_str());i++) {
858 ASSERT(tok_iter != tokens.end(),
859 "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
860 if (_timestep == 0)
861 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
862 lineStream << word << "\t";
863 tok_iter++;
864 }
865 if (_timestep == 0) {
866 readNeighbors(&lineStream,
867 atoi(it->substr(it->find("=") + 1, 1).c_str()), atomid);
868 }
869 break;
870 case TremoloKey::charge :
871 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
872 if (_timestep == 0) {
873 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
874 atomInfo->set(currentField, word);
875 newAtom->setCharge(boost::lexical_cast<double>(word));
876 }
877 tok_iter++;
878 break;
879 default :
880 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
881 if (_timestep == 0) {
882 LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
883 atomInfo->set(currentField, word);
884 }
885 tok_iter++;
886 break;
887 }
888 }
889 LOG(3, "INFO: Parsed atom " << atomid << ".");
890 if (newmol != NULL)
891 newmol->AddAtom(newAtom);
892}
893
894bool FormatParser< tremolo >::saveAtomsInExttypes(
895 std::ostream &output,
896 const std::vector<const atom*> &atoms,
897 const int id) const
898{
899 bool status = true;
900 // parse the file
901 for (std::vector<const atom *>::const_iterator iter = atoms.begin();
902 iter != atoms.end(); ++iter) {
903 const int atomicid = getLocalId((*iter)->getId());
904 if (atomicid == -1)
905 status = false;
906 output << atomicid << "\t" << id << std::endl;
907 }
908
909 return status;
910}
911
912/**
913 * Reads neighbor information for one atom from the input.
914 *
915 * \param line stream where to read the information from
916 * \param numberOfNeighbors number of neighbors to read
917 * \param atomid world id of the atom the information belongs to
918 */
919void FormatParser< tremolo >::readNeighbors(
920 std::stringstream* line,
921 const int numberOfNeighbors,
922 const int atomId) {
923 int neighborId = 0;
924 for (int i = 0; i < numberOfNeighbors; i++) {
925 *line >> neighborId;
926 // 0 is used to fill empty neighbor positions in the tremolo file.
927 if (neighborId > 0) {
928 LOG(4, "INFO: Atom with global id " << atomId
929 << " has neighbour with serial " << neighborId);
930 additionalAtomData[atomId].neighbors.push_back(neighborId);
931 }
932 }
933}
934
935/**
936 * Checks whether the provided name is within \a fields.
937 *
938 * \param fields which usedFields to use
939 * \param fieldName name to check
940 * \return true if the field name is used
941 */
942bool FormatParser< tremolo >::isUsedField(
943 const usedFields_t &fields,
944 const std::string &fieldName) const
945{
946 bool fieldNameExists = false;
947 for (usedFields_t::const_iterator usedField = fields.begin();
948 usedField != fields.end(); usedField++) {
949 if (usedField->substr(0, usedField->find("=")) == fieldName)
950 fieldNameExists = true;
951 }
952
953 return fieldNameExists;
954}
955
956
957/**
958 * Adds the collected neighbor information to the atoms in the world. The atoms
959 * are found by their current ID and mapped to the corresponding atoms with the
960 * Id found in the parsed file.
961 *
962 * @param atoms vector with all newly added (global) atomic ids
963 */
964void FormatParser< tremolo >::processNeighborInformation(
965 const std::vector<atomId_t> &atoms) {
966 if (!isUsedField(usedFields_load, "neighbors")) {
967 return;
968 }
969
970 for (std::vector<atomId_t>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
971 ASSERT(additionalAtomData.count(*iter) != 0,
972 "FormatParser< tremolo >::processNeighborInformation() - global id "
973 +toString(*iter)+" unknown in additionalAtomData.");
974 TremoloAtomInfoContainer &currentInfo = additionalAtomData[*iter];
975 ASSERT (!currentInfo.neighbors_processed,
976 "FormatParser< tremolo >::processNeighborInformation() - neighbors of new atom "
977 +toString(*iter)+" are already processed.");
978 for(std::vector<int>::const_iterator neighbor = currentInfo.neighbors.begin();
979 neighbor != currentInfo.neighbors.end(); neighbor++
980 ) {
981 LOG(3, "INFO: Creating bond between ("
982 << *iter
983 << ") and ("
984 << getGlobalId(*neighbor) << "|" << *neighbor << ")");
985 ASSERT(getGlobalId(*neighbor) != -1,
986 "FormatParser< tremolo >::processNeighborInformation() - global id to local id "
987 +toString(*neighbor)+" is unknown.");
988 World::getInstance().getAtom(AtomById(*iter))
989 ->addBond(WorldTime::getTime(), World::getInstance().getAtom(AtomById(getGlobalId(*neighbor))));
990 }
991 currentInfo.neighbors_processed = true;
992 }
993}
994
995/**
996 * Replaces atom IDs read from the file by the corresponding world IDs. All IDs
997 * IDs of the input string will be replaced; expected separating characters are
998 * "-" and ",".
999 *
1000 * \param string in which atom IDs should be adapted
1001 * \param idgetter function pointer to change the id
1002 *
1003 * \return input string with modified atom IDs
1004 */
1005std::string FormatParser< tremolo >::adaptIdDependentDataString(
1006 const std::string &data,
1007 const boost::function<int (const int)> &idgetter
1008 ) {
1009 // there might be no IDs
1010 if (data == "-") {
1011 return "-";
1012 }
1013
1014 char separator;
1015 int id;
1016 std::stringstream line, result;
1017 line << data;
1018
1019 line >> id;
1020 result << idgetter(id);
1021 while (line.good()) {
1022 line >> separator >> id;
1023 result << separator << idgetter(id);
1024 }
1025
1026 return result.str();
1027}
1028
1029/** Corrects the atom IDs in each imprData entry to the corresponding world IDs
1030 * as they might differ from the originally read IDs.
1031 *
1032 * \param atoms currently parsed in atoms
1033 */
1034void FormatParser< tremolo >::adaptImprData(const std::vector<atomId_t> &atoms) {
1035 if (!isUsedField(usedFields_load, "imprData")) {
1036 return;
1037 }
1038
1039 for (std::vector<atomId_t>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
1040 ASSERT(additionalAtomData.count(*iter) != 0,
1041 "FormatParser< tremolo >::processNeighborInformation() - global id "
1042 +toString(*iter)+" unknown in additionalAtomData.");
1043 TremoloAtomInfoContainer &currentInfo = additionalAtomData[*iter];
1044 currentInfo.imprData = adaptIdDependentDataString(currentInfo.imprData, idglobalizer);
1045 }
1046}
1047
1048/** Corrects the atom IDs in each torsion entry to the corresponding world IDs
1049 * as they might differ from the originally read IDs.
1050 *
1051 * \param atoms currently parsed in atoms
1052 */
1053void FormatParser< tremolo >::adaptTorsion(const std::vector<atomId_t> &atoms) {
1054 if (!isUsedField(usedFields_load, "torsion")) {
1055 return;
1056 }
1057
1058 for (std::vector<atomId_t>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
1059 ASSERT(additionalAtomData.count(*iter) != 0,
1060 "FormatParser< tremolo >::processNeighborInformation() - global id "
1061 +toString(*iter)+" unknown in additionalAtomData.");
1062 TremoloAtomInfoContainer &currentInfo = additionalAtomData[*iter];
1063 currentInfo.torsion = adaptIdDependentDataString(currentInfo.torsion, idglobalizer);
1064 }
1065}
1066
Note: See TracBrowser for help on using the repository browser.