source: src/Parser/TremoloParser.cpp@ a86cda

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

FIX: TremoloParser tripped over following tab in ATOMDATA line.

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