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
RevLine 
[bcf653]1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
[0aa122]4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
[bcf653]5 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
6 */
7
[6bc51d]8/*
9 * TremoloParser.cpp
10 *
11 * Created on: Mar 2, 2010
12 * Author: metzler
13 */
14
[bf3817]15// include config.h
16#ifdef HAVE_CONFIG_H
17#include <config.h>
18#endif
19
[ad011c]20#include "CodePatterns/MemDebug.hpp"
[112b09]21
[ad011c]22#include "CodePatterns/Assert.hpp"
23#include "CodePatterns/Log.hpp"
[4d4d33]24#include "CodePatterns/toString.hpp"
[ad011c]25#include "CodePatterns/Verbose.hpp"
[42127c]26
[9131f3]27#include "TremoloParser.hpp"
[42127c]28
[6f0841]29#include "Atom/atom.hpp"
[129204]30#include "Bond/bond.hpp"
[ccb487]31#include "Box.hpp"
[42127c]32#include "Descriptors/AtomIdDescriptor.hpp"
[3bdb6d]33#include "Element/element.hpp"
34#include "Element/periodentafel.hpp"
[ccb487]35#include "LinearAlgebra/RealSpaceMatrix.hpp"
[42127c]36#include "molecule.hpp"
37#include "MoleculeListClass.hpp"
38#include "World.hpp"
39#include "WorldTime.hpp"
40
[9131f3]41
[9f8b01]42#include <algorithm>
[2034f3]43#include <boost/lexical_cast.hpp>
[72d108]44#include <boost/tokenizer.hpp>
[74a444]45#include <iostream>
46#include <iomanip>
[8bf9c6]47#include <map>
48#include <sstream>
49#include <vector>
[9131f3]50
[5a667d]51#include <boost/assign/list_of.hpp> // for 'map_list_of()'
52#include <boost/assert.hpp>
53
[765f16]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
[5a667d]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
[9131f3]84/**
85 * Constructor.
86 */
[765f16]87FormatParser< tremolo >::FormatParser() :
88 FormatParser_common(NULL)
89{
[4d4d33]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();
[9131f3]97}
98
[5a667d]99
[9131f3]100/**
101 * Destructor.
102 */
[765f16]103FormatParser< tremolo >::~FormatParser()
104{
[23fd43]105 usedFields_save.clear();
[b8d4a3]106 additionalAtomData.clear();
107}
108
109/**
110 * Loads atoms from a tremolo-formatted file.
111 *
112 * \param tremolo file
113 */
[765f16]114void FormatParser< tremolo >::load(istream* file) {
[8bf9c6]115 std::string line;
116 std::string::size_type location;
[b8d4a3]117
[c0e28c]118 // reset the id maps
119 resetIdAssociations();
120
[dc1d9e]121 molecule *newmol = World::getInstance().createMolecule();
[bd2390]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);
[b8d4a3]125 while (file->good()) {
126 std::getline(*file, line, '\n');
[23fd43]127 // we only parse in the first ATOMDATA line
128 if (usedFields_load.empty()) {
[b8d4a3]129 location = line.find("ATOMDATA", 0);
130 if (location != string::npos) {
[23fd43]131 parseAtomDataKeysLine(line, location + 8, usedFields_load);
[b8d4a3]132 }
133 }
134 if (line.length() > 0 && line.at(0) != '#') {
[dc1d9e]135 readAtomDataLine(line, newmol);
[b8d4a3]136 }
137 }
[23fd43]138 LOG(3, "DEBUG: Local usedFields is: " << usedFields_load);
[2e352f]139
[9f8b01]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);
[b8d4a3]145 adaptImprData();
146 adaptTorsion();
[23fd43]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();
[b8d4a3]151}
152
153/**
[73916f]154 * Saves the \a atoms into as a tremolo file.
[b8d4a3]155 *
156 * \param file where to save the state
[73916f]157 * \param atoms atoms to store
[b8d4a3]158 */
[23fd43]159void FormatParser< tremolo >::save(std::ostream* file, const std::vector<atom *> &AtomList) {
[47d041]160 LOG(0, "Saving changes to tremolo.");
[e97a44]161
[23fd43]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);
[b8d4a3]182
[23fd43]183 // store particles
184 for (std::vector<atom*>::const_iterator atomIt = AtomList.begin();
185 atomIt != AtomList.end(); ++atomIt)
186 saveLine(file, *atomIt);
187}
[acd638]188
[5b0581]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
[23fd43]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 */
[27cfde]245void FormatParser< tremolo >::makeUsedFieldsUnique(usedFields_t &fields) const
[23fd43]246{
247 // std::unique only removes if predecessor is equal, not over whole range, hence do it manually
[27cfde]248 usedFields_t temp_fields(fields);
[5b0581]249 usedFieldsSpecialOrderer SpecialOrderer;
250 usedFieldsWeakComparator WeakComparator;
251 std::sort(temp_fields.begin(), temp_fields.end(), SpecialOrderer);
[23fd43]252 usedFields_t::iterator it =
[5b0581]253 std::unique(temp_fields.begin(), temp_fields.end(), WeakComparator);
[23fd43]254 temp_fields.erase(it, temp_fields.end());
[27cfde]255 usedFields_t usedfields(fields);
256 fields.clear();
257 fields.reserve(temp_fields.size());
[23fd43]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()) {
[27cfde]264 fields.push_back(*iter);
[23fd43]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{
[812155]280 resetIdAssociations();
281 atomId_t lastid = 0;
[23fd43]282 for (std::vector<atom*>::const_iterator atomIt = AtomList.begin();
283 atomIt != AtomList.end(); ++atomIt)
[812155]284 associateLocaltoGlobalId(++lastid, (*atomIt)->getId());
[23fd43]285}
[812155]286
[23fd43]287/** Store Atomdata line to \a file.
288 *
289 * @param file output stream
290 */
291void FormatParser< tremolo >::save_AtomDataLine(std::ostream* file) const
292{
[b8d4a3]293 *file << "# ATOMDATA";
[23fd43]294 for (usedFields_t::const_iterator it=usedFields_save.begin();
295 it != usedFields_save.end(); ++it)
[b8d4a3]296 *file << "\t" << *it;
[23fd43]297 *file << std::endl;
298}
[ccb487]299
[23fd43]300/** Store Box info to \a file
301 *
302 * @param file output stream
303 */
304void FormatParser< tremolo >::save_BoxLine(std::ostream* file) const
305{
[ccb487]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;
[b8d4a3]312}
313
[6bc86c]314/** Add default info, when new atom is added to World.
315 *
316 * @param id of atom
317 */
[765f16]318void FormatParser< tremolo >::AtomInserted(atomId_t id)
[6bc86c]319{
[8bf9c6]320 std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
[6bc86c]321 ASSERT(iter == additionalAtomData.end(),
[765f16]322 "FormatParser< tremolo >::AtomInserted() - additionalAtomData already present for newly added atom "
[6bc86c]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 */
[765f16]332void FormatParser< tremolo >::AtomRemoved(atomId_t id)
[6bc86c]333{
[8bf9c6]334 std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
[6bc86c]335 // as we do not insert AtomData on AtomInserted, we cannot be assured of its presence
336// ASSERT(iter != additionalAtomData.end(),
[765f16]337// "FormatParser< tremolo >::AtomRemoved() - additionalAtomData is not present for atom "
[6bc86c]338// +toString(id)+" to remove.");
339 if (iter != additionalAtomData.end())
340 additionalAtomData.erase(iter);
341}
342
[b8d4a3]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 */
[23fd43]349void FormatParser< tremolo >::saveLine(std::ostream* file, const atom* currentAtom)
350{
[b8d4a3]351 TremoloKey::atomDataKey currentField;
352
[47d041]353 LOG(4, "INFO: Saving atom " << *currentAtom << ", its father id is " << currentAtom->GetTrueFather()->getId());
[4d4d33]354
[23fd43]355 for (usedFields_t::iterator it = usedFields_save.begin(); it != usedFields_save.end(); it++) {
[b8d4a3]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
[47d041]360 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getPosition());
[d74077]361 *file << currentAtom->at(0) << "\t";
362 *file << currentAtom->at(1) << "\t";
363 *file << currentAtom->at(2) << "\t";
[b8d4a3]364 break;
365 case TremoloKey::u :
366 // for the moment, assume there are always three dimensions
[47d041]367 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicVelocity());
[bce72c]368 *file << currentAtom->getAtomicVelocity()[0] << "\t";
369 *file << currentAtom->getAtomicVelocity()[1] << "\t";
370 *file << currentAtom->getAtomicVelocity()[2] << "\t";
[b8d4a3]371 break;
[305e7e]372 case TremoloKey::type :
[acd638]373 if (additionalAtomData.count(currentAtom->getId())) {
374 if (additionalAtomData[currentAtom->getId()].get(currentField) != "-") {
[47d041]375 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
[acd638]376 *file << additionalAtomData[currentAtom->getId()].get(currentField) << "\t";
377 } else {
[47d041]378 LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value: " << currentAtom->getType()->getSymbol());
[acd638]379 *file << currentAtom->getType()->getSymbol() << "\t";
380 }
381 } else if (additionalAtomData.count(currentAtom->GetTrueFather()->getId())) {
382 if (additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField) != "-") {
[47d041]383 LOG(3, "Writing for type " << knownKeyNames[currentField] << " stuff from father: " << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField));
[acd638]384 *file << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField) << "\t";
385 } else {
[47d041]386 LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value from father: " << currentAtom->GetTrueFather()->getType()->getSymbol());
[acd638]387 *file << currentAtom->GetTrueFather()->getType()->getSymbol() << "\t";
388 }
[4d4d33]389 } else {
[47d041]390 LOG(3, "Writing for type " << knownKeyNames[currentField] << " its default value: " << currentAtom->getType()->getSymbol());
[acd638]391 *file << currentAtom->getType()->getSymbol() << "\t";
[4d4d33]392 }
[b8d4a3]393 break;
394 case TremoloKey::Id :
[47d041]395 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getId()+1);
[812155]396 *file << getLocalId(currentAtom->getId()) << "\t";
[b8d4a3]397 break;
398 case TremoloKey::neighbors :
[47d041]399 LOG(3, "Writing type " << knownKeyNames[currentField]);
[b8d4a3]400 writeNeighbors(file, atoi(it->substr(it->find("=") + 1, 1).c_str()), currentAtom);
401 break;
[74a444]402 case TremoloKey::resSeq :
[4d4d33]403 if (additionalAtomData.count(currentAtom->getId())) {
[47d041]404 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
[74a444]405 *file << additionalAtomData[currentAtom->getId()].get(currentField);
406 } else if (currentAtom->getMolecule() != NULL) {
[47d041]407 LOG(3, "Writing for type " << knownKeyNames[currentField] << " its own id: " << currentAtom->getMolecule()->getId()+1);
[74a444]408 *file << setw(4) << currentAtom->getMolecule()->getId()+1;
409 } else {
[47d041]410 LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value: " << defaultAdditionalData.get(currentField));
[74a444]411 *file << defaultAdditionalData.get(currentField);
412 }
413 *file << "\t";
[4d4d33]414 break;
[2034f3]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;
[b8d4a3]433 default :
[4d4d33]434 if (additionalAtomData.count(currentAtom->getId())) {
[47d041]435 LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
[74a444]436 *file << additionalAtomData[currentAtom->getId()].get(currentField);
[4d4d33]437 } else if (additionalAtomData.count(currentAtom->GetTrueFather()->getId())) {
[47d041]438 LOG(3, "Writing for type " << knownKeyNames[currentField] << " stuff from father: " << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField));
[74a444]439 *file << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField);
440 } else {
[47d041]441 LOG(3, "Writing for type " << knownKeyNames[currentField] << " the default: " << defaultAdditionalData.get(currentField));
[74a444]442 *file << defaultAdditionalData.get(currentField);
443 }
[b8d4a3]444 *file << "\t";
445 break;
446 }
447 }
448
[23fd43]449 *file << std::endl;
[b8d4a3]450}
451
452/**
453 * Writes the neighbor information of one atom to the provided stream.
454 *
[9d83b6]455 * Note that ListOfBonds of WorldTime::CurrentTime is used.
456 *
[b8d4a3]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 */
[23fd43]461void FormatParser< tremolo >::writeNeighbors(std::ostream* file, const int numberOfNeighbors, const atom* currentAtom) {
[9d83b6]462 const BondList& ListOfBonds = currentAtom->getListOfBonds();
[ca2cfa]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)
[812155]468 sortedBonds.insert(getLocalId((*iter)->GetOtherAtom(currentAtom)->getId()));
[ca2cfa]469 // print indices
470 sortedIndices::const_iterator currentBond = sortedBonds.begin();
[b8d4a3]471 for (int i = 0; i < numberOfNeighbors; i++) {
[812155]472 *file << (currentBond != sortedBonds.end() ? (*currentBond) : 0) << "\t";
[ca2cfa]473 if (currentBond != sortedBonds.end())
[0bbfa1]474 ++currentBond;
[b8d4a3]475 }
[9131f3]476}
477
478/**
[23fd43]479 * Stores keys from the ATOMDATA line in \a fields.
[9131f3]480 *
481 * \param line to parse the keys from
[23fd43]482 * \param offset with which offset the keys begin within the line
483 * \param fields which usedFields to use
[9131f3]484 */
[23fd43]485void FormatParser< tremolo >::parseAtomDataKeysLine(
486 const std::string &line,
487 const int offset,
488 usedFields_t &fields) {
[8bf9c6]489 std::string keyword;
490 std::stringstream lineStream;
[9131f3]491
492 lineStream << line.substr(offset);
493 while (lineStream.good()) {
494 lineStream >> keyword;
[b8d4a3]495 if (knownKeys[keyword.substr(0, keyword.find("="))] == TremoloKey::noKey) {
[ecb799]496 // TODO: throw exception about unknown key
[5a667d]497 cout << "Unknown key: " << keyword.substr(0, keyword.find("=")) << " is not part of the tremolo format specification." << endl;
498 throw IllegalParserKeyException();
[4415da]499 break;
500 }
[23fd43]501 fields.push_back(keyword);
[9131f3]502 }
[23fd43]503 //LOG(1, "INFO: " << fields);
[9131f3]504}
505
[5a667d]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
[81c980b]533/** Sets the properties per atom to print to .data file by parsing line from
534 * \a atomdata_string.
535 *
[23fd43]536 * We just call \sa FormatParser< tremolo >::parseAtomDataKeysLine(), however
537 * we clear FormatParser< tremolo >::usedFields_save.
[81c980b]538 *
539 * @param atomdata_string line to parse with space-separated values
540 */
[765f16]541void FormatParser< tremolo >::setAtomData(const std::string &atomdata_string)
[81c980b]542{
[23fd43]543 usedFields_save.clear();
544 parseAtomDataKeysLine(atomdata_string, 0, usedFields_save);
[81c980b]545}
546
547
[9131f3]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
[dc1d9e]553 * \param *newmol molecule to add atom to
[9131f3]554 */
[23fd43]555void FormatParser< tremolo >::readAtomDataLine(const std::string &line, molecule *newmol = NULL) {
[8bf9c6]556 std::stringstream lineStream;
[4415da]557 atom* newAtom = World::getInstance().createAtom();
[89a31d]558 const atomId_t atomid = newAtom->getId();
559 additionalAtomData[atomid] = TremoloAtomInfoContainer(); // fill with default values
560 TremoloAtomInfoContainer *atomInfo = &additionalAtomData[atomid];
[b8d4a3]561 TremoloKey::atomDataKey currentField;
[72d108]562 ConvertTo<double> toDouble;
563 ConvertTo<int> toInt;
[056e70]564 Vector tempVector;
[72d108]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(),
[765f16]572 "FormatParser< tremolo >::readAtomDataLine - empty string, need at least ' '!");
[72d108]573 tokenizer::iterator tok_iter = tokens.begin();
574 // then associate each token to each file
[23fd43]575 for (usedFields_t::const_iterator it = usedFields_load.begin(); it < usedFields_load.end(); it++) {
[72d108]576 const std::string keyName = it->substr(0, it->find("="));
577 currentField = knownKeys[keyName];
[8bf9c6]578 const std::string word = *tok_iter;
[47d041]579 LOG(4, "INFO: Parsing key " << keyName << " with remaining data " << word);
[4415da]580 switch (currentField) {
[b8d4a3]581 case TremoloKey::x :
[4415da]582 // for the moment, assume there are always three dimensions
[d74077]583 for (int i=0;i<NDIM;i++) {
[765f16]584 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for x["+toString(i)+"]!");
[47d041]585 LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
[72d108]586 newAtom->set(i, toDouble(*tok_iter));
587 tok_iter++;
[d74077]588 }
[4415da]589 break;
[b8d4a3]590 case TremoloKey::u :
[4415da]591 // for the moment, assume there are always three dimensions
[72d108]592 for (int i=0;i<NDIM;i++) {
[765f16]593 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for u["+toString(i)+"]!");
[47d041]594 LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
[056e70]595 tempVector[i] = toDouble(*tok_iter);
[72d108]596 tok_iter++;
597 }
[056e70]598 newAtom->setAtomicVelocity(tempVector);
[4415da]599 break;
[305e7e]600 case TremoloKey::type :
[4d4d33]601 {
[765f16]602 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
[47d041]603 LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
[a275b3]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 }
[4d4d33]613 // put type name into container for later use
614 atomInfo->set(currentField, *tok_iter);
[47d041]615 LOG(4, "INFO: Parsing element " << (*tok_iter) << " as " << element << " according to KnownTypes.");
[72d108]616 tok_iter++;
[4d4d33]617 newAtom->setType(World::getInstance().getPeriode()->FindElement(element));
[b8d4a3]618 ASSERT(newAtom->getType(), "Type was not set for this atom");
[4415da]619 break;
[4d4d33]620 }
[b8d4a3]621 case TremoloKey::Id :
[765f16]622 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
[47d041]623 LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
[c0e28c]624 associateLocaltoGlobalId(toInt(*tok_iter), atomid);
[72d108]625 tok_iter++;
[4415da]626 break;
[b8d4a3]627 case TremoloKey::neighbors :
[72d108]628 for (int i=0;i<atoi(it->substr(it->find("=") + 1, 1).c_str());i++) {
[765f16]629 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
[47d041]630 LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
[72d108]631 lineStream << *tok_iter << "\t";
632 tok_iter++;
633 }
[b8d4a3]634 readNeighbors(&lineStream,
[89a31d]635 atoi(it->substr(it->find("=") + 1, 1).c_str()), atomid);
[9131f3]636 break;
[2034f3]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;
[9131f3]644 default :
[765f16]645 ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
[47d041]646 LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
[72d108]647 atomInfo->set(currentField, *tok_iter);
648 tok_iter++;
[9131f3]649 break;
650 }
651 }
[89a31d]652 LOG(3, "INFO: Parsed atom " << atomid << ".");
653 if (newmol != NULL)
[dc1d9e]654 newmol->AddAtom(newAtom);
[6bc51d]655}
[9131f3]656
[531f27]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
[b8d4a3]672/**
673 * Reads neighbor information for one atom from the input.
674 *
[0bbfa1]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
[b8d4a3]678 */
[23fd43]679void FormatParser< tremolo >::readNeighbors(std::stringstream* line, const int numberOfNeighbors, const int atomId) {
[b8d4a3]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) {
[47d041]685 LOG(4, "INFO: Atom with global id " << atomId
686 << " has neighbour with serial " << neighborId);
[b8d4a3]687 additionalAtomData[atomId].neighbors.push_back(neighborId);
688 }
689 }
690}
[9131f3]691
692/**
[23fd43]693 * Checks whether the provided name is within \a fields.
[b8d4a3]694 *
[23fd43]695 * \param fields which usedFields to use
696 * \param fieldName name to check
[b8d4a3]697 * \return true if the field name is used
[9131f3]698 */
[23fd43]699bool FormatParser< tremolo >::isUsedField(const usedFields_t &fields, const std::string &fieldName) const
700{
[b8d4a3]701 bool fieldNameExists = false;
[23fd43]702 for (usedFields_t::const_iterator usedField = fields.begin();
703 usedField != fields.end(); usedField++) {
[b8d4a3]704 if (usedField->substr(0, usedField->find("=")) == fieldName)
705 fieldNameExists = true;
706 }
[9131f3]707
[b8d4a3]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.
[9f8b01]716 *
717 * @param atoms vector with all newly added (global) atomic ids
[b8d4a3]718 */
[9f8b01]719void FormatParser< tremolo >::processNeighborInformation(const std::vector<atomId_t> &atoms) {
[23fd43]720 if (!isUsedField(usedFields_load, "neighbors")) {
[b8d4a3]721 return;
722 }
723
[9f8b01]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))));
[9131f3]744 }
[9f8b01]745 currentInfo.neighbors_processed = true;
[9131f3]746 }
[6bc51d]747}
748
[9131f3]749/**
[b8d4a3]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 ",".
[9131f3]753 *
[b8d4a3]754 * \param string in which atom IDs should be adapted
755 *
756 * \return input string with modified atom IDs
[9131f3]757 */
[955b91]758std::string FormatParser< tremolo >::adaptIdDependentDataString(std::string data) {
[b8d4a3]759 // there might be no IDs
760 if (data == "-") {
761 return "-";
762 }
763
764 char separator;
765 int id;
[8bf9c6]766 std::stringstream line, result;
[b8d4a3]767 line << data;
768
769 line >> id;
[c0e28c]770 result << getGlobalId(id);
[b8d4a3]771 while (line.good()) {
772 line >> separator >> id;
[c0e28c]773 result << separator << getGlobalId(id);
[b8d4a3]774 }
775
776 return result.str();
[6bc51d]777}
[b8d4a3]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 */
[765f16]783void FormatParser< tremolo >::adaptImprData() {
[23fd43]784 if (!isUsedField(usedFields_load, "imprData")) {
[b8d4a3]785 return;
786 }
787
[8bf9c6]788 for(std::map<const atomId_t, TremoloAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
[b8d4a3]789 currentInfo != additionalAtomData.end(); currentInfo++
790 ) {
791 currentInfo->second.imprData = adaptIdDependentDataString(currentInfo->second.imprData);
792 }
[6bc51d]793}
[4415da]794
[b8d4a3]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 */
[765f16]799void FormatParser< tremolo >::adaptTorsion() {
[23fd43]800 if (!isUsedField(usedFields_load, "torsion")) {
[b8d4a3]801 return;
802 }
803
[8bf9c6]804 for(std::map<const atomId_t, TremoloAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
[b8d4a3]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.