source: src/Parser/TremoloParser.cpp@ 5b0581

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

FIX: TremoloParser did not make usedFields unique when neighbors=4 and neighbors=2 was given.

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