source: src/Parser/TremoloParser.cpp@ 05e2ed

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 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 05e2ed was 05e2ed, checked in by Frederik Heber <heber@…>, 13 years ago

FIX: TremoloParser did not load/save imprData or torsion correctly.

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