source: src/Parser/TremoloParser.cpp@ df5b8c

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

Converted FormatParser::save() to using vector of const atom ptrs.

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