source: src/Parser/TremoloParser.cpp@ f2d5ce

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * TremoloParser.cpp
 *
 *  Created on: Mar 2, 2010
 *      Author: metzler
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

//#include "CodePatterns/MemDebug.hpp"

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/toString.hpp"
#include "CodePatterns/Verbose.hpp"

#include "TremoloParser.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "Parser/Exceptions.hpp"
#include "Potentials/Particles/ParticleRegistry.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "Potentials/Particles/ParticleRegistry.hpp"
#include "molecule.hpp"
#include "World.hpp"
#include "WorldTime.hpp"


#include <algorithm>
#include <boost/bind.hpp>
#include <boost/function.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/tokenizer.hpp>
#include <iostream>
#include <iomanip>
#include <map>
#include <sstream>
#include <string>
#include <vector>

#include <boost/assign/list_of.hpp> // for 'map_list_of()'
#include <boost/assert.hpp>

// declare specialized static variables
const std::string FormatParserTrait<tremolo>::name = "tremolo";
const std::string FormatParserTrait<tremolo>::suffix = "data";
const ParserTypes FormatParserTrait<tremolo>::type = tremolo;

// static instances
std::map<std::string, TremoloKey::atomDataKey> FormatParser<tremolo>::knownKeys =
    boost::assign::map_list_of("x",TremoloKey::x)
    ("u",TremoloKey::u)
    ("F",TremoloKey::F)
    ("stress",TremoloKey::stress)
    ("Id",TremoloKey::Id)
    ("neighbors",TremoloKey::neighbors)
    ("imprData",TremoloKey::imprData)
    ("GroupMeasureTypeNo",TremoloKey::GroupMeasureTypeNo)
    ("type",TremoloKey::type)
    ("extType",TremoloKey::extType)
    ("name",TremoloKey::name)
    ("resName",TremoloKey::resName)
    ("chainID",TremoloKey::chainID)
    ("resSeq",TremoloKey::resSeq)
    ("occupancy",TremoloKey::occupancy)
    ("tempFactor",TremoloKey::tempFactor)
    ("segID",TremoloKey::segID)
    ("Charge",TremoloKey::Charge)
    ("charge",TremoloKey::charge)
    ("GrpTypeNo",TremoloKey::GrpTypeNo)
    ("torsion",TremoloKey::torsion)
    (" ",TremoloKey::noKey); // with this we can detect invalid keys

/**
 * Constructor.
 */
FormatParser< tremolo >::FormatParser()  :
  FormatParser_common(NULL),
  idglobalizer(boost::bind(&FormatParser< tremolo >::getGlobalId, this, _1)),
  idlocalizer(boost::bind(&FormatParser< tremolo >::getLocalId, this, _1))
{
  // invert knownKeys for debug output
  for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin(); iter != knownKeys.end(); ++iter)
    knownKeyNames.insert( make_pair( iter->second, iter->first) );

  additionalAtomData.clear();
}


/**
 * Destructor.
 */
FormatParser< tremolo >::~FormatParser() 
{
  usedFields_save.clear();
  additionalAtomData.clear();
}

/**
 * Loads atoms from a tremolo-formatted file.
 *
 * \param tremolo file
 */
void FormatParser< tremolo >::load(istream* file) {
  std::string line;
  std::string::size_type location;

  // reset the id maps
  resetIdAssociations();

  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  usedFields_t usedFields_temp;
  size_t timestep = 0;
  atom *addedatom = NULL;
  std::list<atom *> addedatoms;
  std::list<atom *>::iterator atomiter = addedatoms.begin();
  while (file->good()) {
    std::getline(*file, line, '\n');
    // we only parse in the first ATOMDATA line
    location = line.find("ATOMDATA", 0);
    if (location != string::npos) {
      parseAtomDataKeysLine(line, location + 8, usedFields_temp);
      if (usedFields_load.empty()) {
        // first ATOMDATA: use this line
        usedFields_load.insert(usedFields_load.end(), usedFields_temp.begin(), usedFields_temp.end());
        LOG(3, "DEBUG: Local usedFields is: " << usedFields_load);
      } else {
        // following ATOMDATA: check against present
        LOG(3, "DEBUG: Parsed check usedFields is: " << usedFields_temp);
        if (usedFields_load != usedFields_temp) {
          ELOG(1, "File contains multiple time steps with differing ATOMDATA line, parsing only first time step.");
          break;
        } else {
          // update time step
          ++timestep;
          atomiter = addedatoms.begin();
        }
      }
      usedFields_temp.clear();
    } else {
      if (line.length() > 0 && line.at(0) != '#') {
        if (timestep != 0) {
          ASSERT( atomiter != addedatoms.end(),
              "FormatParser< tremolo >::load() - timesteps contains more atoms than first one.");
          addedatom = *atomiter++;
        }
        readAtomDataLine(line, newmol, timestep, addedatom);
        if (timestep == 0)
          addedatoms.push_back(addedatom);
      }
    }
  }
  ASSERT( addedatoms.size() == newmol->size(),
      "FormatParser< tremolo >::load() - number of atoms in mol and addedatoms differ.");

  // refresh atom::nr and atom::name
  std::vector<atomId_t> atoms(newmol->getAtomCount());
  std::transform(newmol->begin(), newmol->end(), atoms.begin(),
      boost::bind(&atom::getId, _1));
  processNeighborInformation(atoms);
  adaptImprData(atoms);
  adaptTorsion(atoms);

  // append usedFields to global usedFields, is made unique on save, clear after use
  usedFields_save.insert(usedFields_save.end(), usedFields_load.begin(), usedFields_load.end());
  usedFields_load.clear();
}

/**
 * Saves the \a atoms into as a tremolo file.
 *
 * \param file where to save the state
 * \param atoms atoms to store
 */
void FormatParser< tremolo >::save(
    std::ostream* file,
    const std::vector<const atom *> &AtomList) {
  LOG(2, "DEBUG: Saving changes to tremolo.");

  // install default usedFields if empty so far
  if (usedFields_save.empty()) {
    // default behavior: use all possible keys on output
    for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin();
        iter != knownKeys.end(); ++iter)
      if (iter->second != TremoloKey::noKey) // don't add noKey
        usedFields_save.push_back(iter->first);
  }
  // make present usedFields_save unique
  makeUsedFieldsUnique(usedFields_save);
  LOG(1, "INFO: Global (with unique entries) usedFields_save is: " << usedFields_save);

  // distribute ids continuously
  distributeContinuousIds(AtomList);

  std::pair<size_t, size_t> minmax_trajectories =
      getMinMaxTrajectories(AtomList);
  LOG(2, "INFO: There are " << minmax_trajectories.second << " steps to save.");

  for (size_t step = 0; (step < minmax_trajectories.second) || (step == 0); ++step) {
    // store atomdata
    save_AtomDataLine(file);

    // store box only on first step
    if (step == 0)
      save_BoxLine(file);

    // store particles
    for (std::vector<const atom*>::const_iterator  atomIt = AtomList.begin();
        atomIt != AtomList.end(); ++atomIt)
      saveLine(file, *atomIt, step);
  }
}

struct usedFieldsWeakComparator
{
  /** Special comparator regards "neighbors=4" and "neighbors=2" as equal
   *
   * \note This one is used for making usedFields unique, i.e. throwing out the "smaller"
   * neighbors.
   */
  bool operator()(const std::string &a, const std::string &b) const
  {
    // only compare up to first equality sign
    return (a.substr(0, a.find_first_of('=')) == b.substr(0, b.find_first_of('=')));
  }
};

struct usedFieldsSpecialOrderer
{
  /** Special string comparator that regards "neighbors=4" < "neighbors=2" as true and
   * the other way round as false.
   *
   * Here, we implement the operator "\a < \b" in a special way to allow the
   * above.
   *
   * \note This one is used for sorting usedFields in preparation for making it unique.
   */
  bool operator()(const std::string &a, const std::string &b) const
  {
    // only compare up to first equality sign
    size_t a_equality = a.find_first_of('=');
    size_t b_equality = b.find_first_of('=');
    // if key before equality is not equal, return whether it is smaller or not
    if (a.substr(0, a_equality) != b.substr(0, b_equality)) {
      return a.substr(0, a_equality) < b.substr(0, b_equality);
    } else { // now we know that the key before equality is the same in either string
      // if one of them has no equality, the one with equality must go before
      if ((a_equality != std::string::npos) && (b_equality == std::string::npos))
        return true;
      if ((a_equality == std::string::npos) && (b_equality != std::string::npos))
        return false;
      // if both don't have equality (and the token before is equal), it is not "<" but "=="
      if ((a_equality == std::string::npos) && (b_equality == std::string::npos))
        return false;
      // if now both have equality sign, the larger value after it, must come first
      return a.substr(a_equality, std::string::npos) > b.substr(b_equality, std::string::npos);
    }
  }
};

/** Helper function to make \given fields unique while preserving the order of first appearance.
 *
 * As std::unique only removes element if equal to predecessor, a vector is only
 * made unique if sorted beforehand. But sorting would destroy order of first
 * appearance, hence we do the sorting on a temporary field and add the unique
 * elements in the order as in \a fields.
 *
 * @param fields usedFields to make unique while preserving order of appearance
 */
void FormatParser< tremolo >::makeUsedFieldsUnique(usedFields_t &fields) const
{
  // std::unique only removes if predecessor is equal, not over whole range, hence do it manually
  usedFields_t temp_fields(fields);
  usedFieldsSpecialOrderer SpecialOrderer;
  usedFieldsWeakComparator WeakComparator;
  std::sort(temp_fields.begin(), temp_fields.end(), SpecialOrderer);
  usedFields_t::iterator it =
      std::unique(temp_fields.begin(), temp_fields.end(), WeakComparator);
  temp_fields.erase(it, temp_fields.end());
  usedFields_t usedfields(fields);
  fields.clear();
  fields.reserve(temp_fields.size());
  // now go through each usedFields entry, check if in temp_fields and remove there on first occurence
  for (usedFields_t::const_iterator iter = usedfields.begin();
      iter != usedfields.end(); ++iter) {
    usedFields_t::iterator uniqueiter =
        std::find(temp_fields.begin(), temp_fields.end(), *iter);
    if (uniqueiter != temp_fields.end()) {
      fields.push_back(*iter);
      // add only once to ATOMDATA
      temp_fields.erase(uniqueiter);
    }
  }
  ASSERT( temp_fields.empty(),
      "FormatParser< tremolo >::save() - still unique entries left in temp_fields after unique?");
}


/** Resets and distributes the indices continuously.
 *
 * \param atoms atoms to store
 */
void FormatParser< tremolo >::distributeContinuousIds(
    const std::vector<const atom *> &AtomList)
{
  resetIdAssociations();
  atomId_t lastid = 0;
  for (std::vector<const atom*>::const_iterator atomIt = AtomList.begin();
      atomIt != AtomList.end(); ++atomIt)
    associateLocaltoGlobalId(++lastid, (*atomIt)->getId());
}

/** Store Atomdata line to \a file.
 *
 * @param file output stream
 */
void FormatParser< tremolo >::save_AtomDataLine(std::ostream* file) const
{
  *file << "# ATOMDATA";
  for (usedFields_t::const_iterator it=usedFields_save.begin();
      it != usedFields_save.end(); ++it)
    *file << "\t" << *it;
  *file << std::endl;
}

/** Store Box info to \a file
 *
 * @param file output stream
 */
void FormatParser< tremolo >::save_BoxLine(std::ostream* file) const
{
  *file << "# Box";
  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
  for (size_t i=0; i<NDIM;++i)
    for (size_t j=0; j<NDIM;++j)
      *file << "\t" << M.at(i,j);
  *file << std::endl;
}

/** Add default info, when new atom is added to World.
 *
 * @param id of atom
 */
void FormatParser< tremolo >::AtomInserted(atomId_t id)
{
  std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
  ASSERT(iter == additionalAtomData.end(),
      "FormatParser< tremolo >::AtomInserted() - additionalAtomData already present for newly added atom "
      +toString(id)+".");
  // don't add entry, as this gives a default resSeq of 0 not the molecule id
  // additionalAtomData.insert( std::make_pair(id, TremoloAtomInfoContainer()) );
}

/** Remove additional AtomData info, when atom has been removed from World.
 *
 * @param id of atom
 */
void FormatParser< tremolo >::AtomRemoved(atomId_t id)
{
  std::map<const atomId_t, TremoloAtomInfoContainer>::iterator iter = additionalAtomData.find(id);
  // as we do not insert AtomData on AtomInserted, we cannot be assured of its presence
//  ASSERT(iter != additionalAtomData.end(),
//      "FormatParser< tremolo >::AtomRemoved() - additionalAtomData is not present for atom "
//      +toString(id)+" to remove.");
  if (iter != additionalAtomData.end())
    additionalAtomData.erase(iter);
}

template <class T>
T NoOp(const atom * const)
{
  return T();
}

template <class T>
void writeEntryFromAdditionalAtomData_ifpresent(
    std::map<const atomId_t, TremoloAtomInfoContainer> &_additionalAtomData,
    std::map<TremoloKey::atomDataKey, std::string> &_knownKeyNames,
    std::ostream* _file,
    const atom * const _currentAtom,
    const TremoloKey::atomDataKey _currentField,
    const typename boost::function<T (const atom * const)> _getter)
{
  if (_additionalAtomData.count(_currentAtom->getId())) {
    LOG(3, "Writing for type " << _knownKeyNames[_currentField] << ": "
        << _additionalAtomData[_currentAtom->getId()].get(_currentField));
    *_file << _additionalAtomData[_currentAtom->getId()].get(_currentField);
  } else if (_additionalAtomData.count(_currentAtom->GetTrueFather()->getId())) {
    LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " stuff from father: "
        << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField));
    *_file << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField);
  } else {
    LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " its default value: "
        << _getter(_currentAtom));
    *_file << _getter(_currentAtom);
  }
  *_file << "\t";
}

template <class T>
void writeAdditionalAtomDataEntry(
    std::map<const atomId_t, TremoloAtomInfoContainer> &_additionalAtomData,
    std::map<TremoloKey::atomDataKey, std::string> &_knownKeyNames,
    TremoloAtomInfoContainer &_defaultAdditionalData,
    std::ostream* _file,
    const atom * const _currentAtom,
    const TremoloKey::atomDataKey _currentField)
{
  if (_additionalAtomData.count(_currentAtom->getId())) {
    LOG(3, "Writing for type " << _knownKeyNames[_currentField] << ": "
        << _additionalAtomData[_currentAtom->getId()].get(_currentField));
    *_file << _additionalAtomData[_currentAtom->getId()].get(_currentField);
  } else if (_additionalAtomData.count(_currentAtom->GetTrueFather()->getId())) {
    LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " stuff from father: "
        << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField));
    *_file << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField);
  } else {
    LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " its default value: "
        << _defaultAdditionalData.get(_currentField));
    *_file << _defaultAdditionalData.get(_currentField);
  }
  *_file << "\t";
}

template <class T>
void writeEntryFromAdditionalAtomData_ifnotempty(
    std::map<const atomId_t, TremoloAtomInfoContainer> &_additionalAtomData,
    std::map<TremoloKey::atomDataKey, std::string> &_knownKeyNames,
    TremoloAtomInfoContainer &_defaultAdditionalData,
    std::ostream* _file,
    const atom * const _currentAtom,
    const TremoloKey::atomDataKey _currentField,
    const typename boost::function<T (const atom * const)> _getter)
{
  if (_additionalAtomData.count(_currentAtom->getId())) {
    if (_additionalAtomData[_currentAtom->getId()].get(_currentField) != "-") {
      LOG(3, "Writing for type " << _knownKeyNames[_currentField] << ": "
          << _additionalAtomData[_currentAtom->getId()].get(_currentField));
      *_file << _additionalAtomData[_currentAtom->getId()].get(_currentField);
    } else {
      LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " default value: "
          << _getter(_currentAtom));
      *_file << _getter(_currentAtom);
    }
  } else if (_additionalAtomData.count(_currentAtom->GetTrueFather()->getId())) {
    if (_additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField) != "-") {
      LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " stuff from father: "
          << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField));
      *_file << _additionalAtomData[_currentAtom->GetTrueFather()->getId()].get(_currentField);
    } else {
      LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " default value from father: "
          << _getter(_currentAtom->GetTrueFather()));
      *_file << _getter(_currentAtom->GetTrueFather());
    }
  } else {
    LOG(3, "Writing for type " << _knownKeyNames[_currentField] << " its default value: "
        << _getter(_currentAtom));
    *_file << _getter(_currentAtom);
  }
  *_file << "\t";
}

/**
 * Writes one line of tremolo-formatted data to the provided stream.
 *
 * \param stream where to write the line to
 * \param reference to the atom of which information should be written
 */
void FormatParser< tremolo >::saveLine(
    std::ostream* file,
    const atom * const currentAtom,
    const size_t _timestep)
{
  TremoloKey::atomDataKey currentField;

  LOG(4, "INFO: Saving atom " << *currentAtom << ", its father id is "
      << currentAtom->GetTrueFather()->getId() << " at time step " << _timestep);

  // note down default precision
  std::streamsize defaultprecision=file->precision();
  std::streamsize highprecision=10;

  for (usedFields_t::iterator it = usedFields_save.begin(); it != usedFields_save.end(); it++) {
    currentField = knownKeys[it->substr(0, it->find("="))];
    switch (currentField) {
      case TremoloKey::x :
        // for the moment, assume there are always three dimensions
        LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getPositionAtStep(_timestep));
        file->precision(highprecision);
        *file << currentAtom->atStep(0, _timestep) << "\t";
        *file << currentAtom->atStep(1, _timestep) << "\t";
        *file << currentAtom->atStep(2, _timestep) << "\t";
        file->precision(defaultprecision);
        break;
      case TremoloKey::u :
        // for the moment, assume there are always three dimensions
        LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicVelocityAtStep(_timestep));
        file->precision(highprecision);
        *file << currentAtom->getAtomicVelocityAtStep(_timestep)[0] << "\t";
        *file << currentAtom->getAtomicVelocityAtStep(_timestep)[1] << "\t";
        *file << currentAtom->getAtomicVelocityAtStep(_timestep)[2] << "\t";
        file->precision(defaultprecision);
        break;
      case TremoloKey::F :
        // for the moment, assume there are always three dimensions
        LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicForceAtStep(_timestep));
        file->precision(highprecision);
        *file << currentAtom->getAtomicForceAtStep(_timestep)[0] << "\t";
        *file << currentAtom->getAtomicForceAtStep(_timestep)[1] << "\t";
        *file << currentAtom->getAtomicForceAtStep(_timestep)[2] << "\t";
        file->precision(defaultprecision);
        break;
      case TremoloKey::type :
        writeEntryFromAdditionalAtomData_ifnotempty<std::string>(
            additionalAtomData,
            knownKeyNames,
            defaultAdditionalData,
            file,
            currentAtom,
            currentField,
            boost::bind(&element::getSymbol, boost::bind(&AtomInfo::getType, _1)));
        break;
      case TremoloKey::Id :
        LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getId()+1);
        *file << getLocalId(currentAtom->getId()) << "\t";
        break;
      case TremoloKey::neighbors :
        LOG(3, "Writing type " << knownKeyNames[currentField]);
        writeNeighbors(file, atoi(it->substr(it->find("=") + 1, 1).c_str()), currentAtom);
        break;
      case TremoloKey::imprData :
      case TremoloKey::torsion :
        LOG(3, "Writing type " << knownKeyNames[currentField]);
        *file << adaptIdDependentDataString(
            additionalAtomData[currentAtom->getId()].get(currentField),
            idlocalizer)
            << "\t";
        break;
      case TremoloKey::resSeq :
        if (additionalAtomData.count(currentAtom->getId())) {
          LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << additionalAtomData[currentAtom->getId()].get(currentField));
          *file << additionalAtomData[currentAtom->getId()].get(currentField);
        } else if (currentAtom->getMolecule() != NULL) {
          LOG(3, "Writing for type " << knownKeyNames[currentField] << " its own id: " << currentAtom->getMolecule()->getId()+1);
          *file << setw(4) << currentAtom->getMolecule()->getId();
        } else {
          LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value: " << defaultAdditionalData.get(currentField));
          *file << defaultAdditionalData.get(currentField);
        }
        *file << "\t";
        break;
      case TremoloKey::charge :
        if (currentAtom->getCharge() == 0.) {
          writeEntryFromAdditionalAtomData_ifpresent<double>(
              additionalAtomData,
              knownKeyNames,
              file,
              currentAtom,
              currentField,
              boost::bind(&AtomInfo::getCharge, _1));
        } else {
          LOG(3, "Writing for type " << knownKeyNames[currentField] << " AtomInfo::charge : " << currentAtom->getCharge());
          *file << currentAtom->getCharge() << "\t";
        }
        break;
      default :
        writeAdditionalAtomDataEntry<std::string>(
            additionalAtomData,
            knownKeyNames,
            defaultAdditionalData,
            file,
            currentAtom,
            currentField);
        break;
    }
  }

  *file << std::endl;
}

/**
 * Writes the neighbor information of one atom to the provided stream.
 *
 * Note that ListOfBonds of WorldTime::CurrentTime is used.
 *
 * \param stream where to write neighbor information to
 * \param number of neighbors
 * \param reference to the atom of which to take the neighbor information
 */
void FormatParser< tremolo >::writeNeighbors(
    std::ostream* file,
    const int numberOfNeighbors,
    const atom * const currentAtom) {
  const BondList& ListOfBonds = currentAtom->getListOfBonds();
  // sort bonded indices
  typedef std::set<atomId_t> sortedIndices;
  sortedIndices sortedBonds;
  for (BondList::const_iterator iter = ListOfBonds.begin();
      iter != ListOfBonds.end(); ++iter)
    sortedBonds.insert(getLocalId((*iter)->GetOtherAtom(currentAtom)->getId()));
  // print indices
  sortedIndices::const_iterator currentBond = sortedBonds.begin();
  for (int i = 0; i < numberOfNeighbors; i++) {
    *file << (currentBond != sortedBonds.end() ? (*currentBond) : 0) << "\t";
    if (currentBond != sortedBonds.end())
      ++currentBond;
  }
}

/**
 * Stores keys from the ATOMDATA line in \a fields.
 *
 * \param line to parse the keys from
 * \param offset with which offset the keys begin within the line
 * \param fields which usedFields to use
 */
void FormatParser< tremolo >::parseAtomDataKeysLine(
    const std::string &line,
    const int offset,
    usedFields_t &fields) {
  std::string keyword;
  std::stringstream lineStream;

  lineStream << line.substr(offset);
  lineStream >> ws;
  while (lineStream.good()) {
    lineStream >> keyword;
    if (knownKeys[keyword.substr(0, keyword.find("="))] == TremoloKey::noKey) {
      // TODO: throw exception about unknown key
      cout << "Unknown key: " << keyword.substr(0, keyword.find("=")) << " is not part of the tremolo format specification." << endl;
      throw IllegalParserKeyException();
      break;
    }
    fields.push_back(keyword);
    lineStream >> ws;
  }
  LOG(2, "INFO: " << fields);
}

/**
 * Tests whether the keys from the ATOMDATA line can be read correctly.
 *
 * \param line to parse the keys from
 */
bool FormatParser< tremolo >::testParseAtomDataKeysLine(
    const std::string &line) {
  std::string keyword;
  std::stringstream lineStream;

  // check string after ATOMDATA
  const std::string AtomData("ATOMDATA");
  const size_t AtomDataOffset = line.find(AtomData, 0);
  if (AtomDataOffset == std::string::npos)
    lineStream << line;
  else
    lineStream << line.substr(AtomDataOffset + AtomData.length());
  while (lineStream.good()) {
    lineStream >> keyword;
    //LOG(2, "DEBUG: Checking key " << keyword.substr(0, keyword.find("=")) << ".");
    if (knownKeys[keyword.substr(0, keyword.find("="))] == TremoloKey::noKey)
      return false;
  }
  //LOG(1, "INFO: " << fields);
  return true;
}

std::string FormatParser< tremolo >::getAtomData() const
{
  std::stringstream output;
  std::for_each(usedFields_save.begin(), usedFields_save.end(),
      output << boost::lambda::_1 << " ");
  const std::string returnstring(output.str());
  return returnstring.substr(0, returnstring.find_last_of(" "));
}

/** Appends the properties per atom to print to .data file by parsing line from
 *  \a atomdata_string.
 *
 *  We just call \sa  FormatParser< tremolo >::parseAtomDataKeysLine().
 *
 * @param atomdata_string line to parse with space-separated values
 */
void FormatParser< tremolo >::setAtomData(const std::string &atomdata_string)
{
  parseAtomDataKeysLine(atomdata_string, 0, usedFields_save);
}

/** Sets the properties per atom to print to .data file by parsing line from
 *  \a atomdata_string.
 *
 *  We just call \sa  FormatParser< tremolo >::parseAtomDataKeysLine(), however
 *  we clear FormatParser< tremolo >::usedFields_save.
 *
 * @param atomdata_string line to parse with space-separated values
 */
void FormatParser< tremolo >::resetAtomData(const std::string &atomdata_string)
{
  usedFields_save.clear();
  parseAtomDataKeysLine(atomdata_string, 0, usedFields_save);
}


/**
 * Reads one data line of a tremolo file and interprets it according to the keys
 * obtained from the ATOMDATA line.
 *
 * \param line to parse as an atom
 * \param *newmol molecule to add atom to
 * \param _timestep time step to parse to
 * \param _addedatom ref to added atom
 */
void FormatParser< tremolo >::readAtomDataLine(
    const std::string &line,
    molecule *newmol,
    const size_t _timestep,
    atom *&_addedatoms)
{
  std::stringstream lineStream;
  atom* newAtom = NULL;
  atomId_t atomid = -1;
  if (_timestep == 0) {
    newAtom = World::getInstance().createAtom();
    _addedatoms = newAtom;
    atomid = newAtom->getId();
    additionalAtomData[atomid] = TremoloAtomInfoContainer(); // fill with default values
  } else {
    newAtom = _addedatoms;
    atomid = newAtom->getId();
  }
  TremoloAtomInfoContainer *atomInfo = &additionalAtomData[atomid];
  TremoloKey::atomDataKey currentField;
  ConvertTo<double> toDouble;
  ConvertTo<int> toInt;
  Vector tempVector;
  const ParticleRegistry& registry = ParticleRegistry::getConstInstance();
  const periodentafel& periode = *World::getInstance().getPeriode();

  // setup tokenizer, splitting up white-spaced entries
  typedef boost::tokenizer<boost::char_separator<char> >
      tokenizer;
  boost::char_separator<char> whitespacesep(" \t");
  tokenizer tokens(line, whitespacesep);
  ASSERT(tokens.begin() != tokens.end(),
      "FormatParser< tremolo >::readAtomDataLine - empty string, need at least ' '!");
  tokenizer::const_iterator tok_iter = tokens.begin();
  // then associate each token to each file
  for (usedFields_t::const_iterator it = usedFields_load.begin(); it != usedFields_load.end(); it++) {
    const std::string keyName = it->substr(0, it->find("="));
    currentField = knownKeys[keyName];
    ASSERT(tok_iter != tokens.end(),
        "FormatParser< tremolo >::readAtomDataLine - too few entries in line '"+line+"'!");
    const std::string &word = *tok_iter;
    LOG(4, "INFO: Parsing key " << keyName << " with remaining data " << word);
    switch (currentField) {
      case TremoloKey::x :
        // for the moment, assume there are always three dimensions
        for (int i=0;i<NDIM;i++) {
          ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for x["+toString(i)+"]!");
          LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
          newAtom->setAtStep(i, _timestep, toDouble(word));
          tok_iter++;
        }
        break;
      case TremoloKey::u :
        // for the moment, assume there are always three dimensions
        for (int i=0;i<NDIM;i++) {
          ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for u["+toString(i)+"]!");
          LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
          tempVector[i] = toDouble(word);
          tok_iter++;
        }
        newAtom->setAtomicVelocityAtStep(_timestep, tempVector);
        break;
      case TremoloKey::F :
        // for the moment, assume there are always three dimensions
        for (int i=0;i<NDIM;i++) {
          ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for F["+toString(i)+"]!");
          LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
          tempVector[i] = toDouble(word);
          tok_iter++;
        }
        newAtom->setAtomicForceAtStep(_timestep, tempVector);
        break;
      case TremoloKey::type :
      {
        ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
        LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
        std::string elementname;
        const element * elem = NULL;
        if (!registry.isPresentByName(word)) {
          std::string lowercase_word(word);
          std::transform(word.begin()+1, word.end(), lowercase_word.begin()+1, ::tolower);
          elem = periode.FindElement(lowercase_word);
          if (elem == NULL) {
            // clean up
            World::getInstance().destroyAtom(newAtom);
            // give an error
            ELOG(0, "TremoloParser: tokens " << word << "/" << lowercase_word
                << " is unknown to neither ParticleRegistry nor Periodentafel.");
            return;
          } else {
            elementname = elem->getSymbol();
          }
        } else {
          const Particle * const p = registry.getByName(word);
          elementname = p->getElement();
          elem = periode.FindElement(elementname);
        }
        // put type name into container for later use
        atomInfo->set(currentField, word);
        LOG(4, "INFO: Parsing element " << (word) << " as " << elementname << " according to KnownTypes.");
        tok_iter++;
        if (newAtom->getType() != elem)
          newAtom->setType(elem);
        ASSERT(newAtom->getType(), "Type was not set for this atom");
        break;
      }
      case TremoloKey::Id :
        ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
        LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
        if (_timestep == 0) {
          associateLocaltoGlobalId(toInt(word), atomid);
        } else {
          // check association is the same
          ASSERT( (atomId_t)getGlobalId(toInt(word)) == atomid,
              "FormatParser< tremolo >::readAtomDataLine() - differing global id "+toString(atomid)
              +" on timestep "+toString(_timestep));
          ASSERT( getLocalId(atomid) == toInt(word),
              "FormatParser< tremolo >::readAtomDataLine() - differing local id "+toString(toInt(word))
              +" on timestep "+toString(_timestep));
        }
        tok_iter++;
        break;
      case TremoloKey::neighbors :
        for (int i=0;i<atoi(it->substr(it->find("=") + 1, 1).c_str());i++) {
          ASSERT(tok_iter != tokens.end(),
              "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
          if (_timestep == 0)
            LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
          lineStream << word << "\t";
          tok_iter++;
        }
        if (_timestep == 0) {
          readNeighbors(&lineStream,
              atoi(it->substr(it->find("=") + 1, 1).c_str()), atomid);
        }
        break;
      case TremoloKey::charge :
        ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
        if (_timestep == 0) {
          LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
          atomInfo->set(currentField, word);
          const double newCharge = boost::lexical_cast<double>(word);
          if (newAtom->getCharge() != newCharge)
            newAtom->setCharge(newCharge);
        }
        tok_iter++;
        break;
      default :
        ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
        if (_timestep == 0) {
          LOG(4, "INFO: Parsing key " << keyName << " with next token " << word << " on step " << _timestep);
          atomInfo->set(currentField, word);
        }
        tok_iter++;
        break;
    }
  }
  LOG(3, "INFO: Parsed atom " << atomid << ".");
  if ((newmol != NULL) && (_timestep == 0))
    newmol->AddAtom(newAtom);
}

bool FormatParser< tremolo >::saveAtomsInExttypes(
    std::ostream &output,
    const std::vector<const atom*> &atoms,
    const int id) const
{
  bool status = true;
  // parse the file
  for (std::vector<const atom *>::const_iterator iter = atoms.begin();
      iter != atoms.end(); ++iter) {
    const int atomicid = getLocalId((*iter)->getId());
    if (atomicid == -1)
      status = false;
    output << atomicid << "\t" << id << std::endl;
  }

  return status;
}

/**
 * Reads neighbor information for one atom from the input.
 *
 * \param line stream where to read the information from
 * \param numberOfNeighbors number of neighbors to read
 * \param atomid world id of the atom the information belongs to
 */
void FormatParser< tremolo >::readNeighbors(
    std::stringstream* line,
    const int numberOfNeighbors,
    const int atomId) {
  int neighborId = 0;
  for (int i = 0; i < numberOfNeighbors; i++) {
    *line >> neighborId;
    // 0 is used to fill empty neighbor positions in the tremolo file.
    if (neighborId > 0) {
      LOG(4, "INFO: Atom with global id " << atomId
          << " has neighbour with serial " << neighborId);
      additionalAtomData[atomId].neighbors.push_back(neighborId);
    }
  }
}

/**
 * Checks whether the provided name is within \a fields.
 *
 * \param fields which usedFields to use
 * \param fieldName name to check
 * \return true if the field name is used
 */
bool FormatParser< tremolo >::isUsedField(
    const usedFields_t &fields,
    const std::string &fieldName) const
{
  bool fieldNameExists = false;
  for (usedFields_t::const_iterator usedField = fields.begin();
      usedField != fields.end(); usedField++) {
    if (usedField->substr(0, usedField->find("=")) == fieldName)
      fieldNameExists = true;
  }

  return fieldNameExists;
}


/**
 * Adds the collected neighbor information to the atoms in the world. The atoms
 * are found by their current ID and mapped to the corresponding atoms with the
 * Id found in the parsed file.
 *
 * @param atoms vector with all newly added (global) atomic ids
 */
void FormatParser< tremolo >::processNeighborInformation(
    const std::vector<atomId_t> &atoms) {
  if (!isUsedField(usedFields_load, "neighbors")) {
    return;
  }

  for (std::vector<atomId_t>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    ASSERT(additionalAtomData.count(*iter) != 0,
        "FormatParser< tremolo >::processNeighborInformation() - global id "
        +toString(*iter)+" unknown in additionalAtomData.");
    TremoloAtomInfoContainer &currentInfo = additionalAtomData[*iter];
    ASSERT (!currentInfo.neighbors_processed,
        "FormatParser< tremolo >::processNeighborInformation() - neighbors of new atom "
        +toString(*iter)+" are already processed.");
    for(std::vector<int>::const_iterator neighbor = currentInfo.neighbors.begin();
      neighbor != currentInfo.neighbors.end(); neighbor++
    ) {
      LOG(3, "INFO: Creating bond between ("
          << *iter
          << ") and ("
          << getGlobalId(*neighbor) << "|" << *neighbor << ")");
      ASSERT(getGlobalId(*neighbor) != -1,
          "FormatParser< tremolo >::processNeighborInformation() - global id to local id "
          +toString(*neighbor)+" is unknown.");
      World::getInstance().getAtom(AtomById(*iter))
          ->addBond(WorldTime::getTime(), World::getInstance().getAtom(AtomById(getGlobalId(*neighbor))));
    }
    currentInfo.neighbors_processed = true;
  }
}

/**
 * Replaces atom IDs read from the file by the corresponding world IDs. All IDs
 * IDs of the input string will be replaced; expected separating characters are
 * "-" and ",".
 *
 * \param string in which atom IDs should be adapted
 * \param idgetter function pointer to change the id
 *
 * \return input string with modified atom IDs
 */
std::string FormatParser< tremolo >::adaptIdDependentDataString(
    const std::string &data,
    const boost::function<int (const int)> &idgetter
    ) {
  // there might be no IDs
  if (data == "-") {
    return "-";
  }

  char separator;
  int id;
  std::stringstream line, result;
  line << data;

  line >> id;
  result << idgetter(id);
  while (line.good()) {
    line >> separator >> id;
    result << separator << idgetter(id);
  }

  return result.str();
}

/** Corrects the atom IDs in each imprData entry to the corresponding world IDs
 * as they might differ from the originally read IDs.
 *
 * \param atoms currently parsed in atoms
 */
void FormatParser< tremolo >::adaptImprData(const std::vector<atomId_t> &atoms) {
  if (!isUsedField(usedFields_load, "imprData")) {
    return;
  }

  for (std::vector<atomId_t>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    ASSERT(additionalAtomData.count(*iter) != 0,
        "FormatParser< tremolo >::processNeighborInformation() - global id "
        +toString(*iter)+" unknown in additionalAtomData.");
    TremoloAtomInfoContainer &currentInfo = additionalAtomData[*iter];
    currentInfo.imprData = adaptIdDependentDataString(currentInfo.imprData, idglobalizer);
  }
}

/** Corrects the atom IDs in each torsion entry to the corresponding world IDs
 * as they might differ from the originally read IDs.
 *
 * \param atoms currently parsed in atoms
 */
void FormatParser< tremolo >::adaptTorsion(const std::vector<atomId_t> &atoms) {
  if (!isUsedField(usedFields_load, "torsion")) {
    return;
  }

  for (std::vector<atomId_t>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    ASSERT(additionalAtomData.count(*iter) != 0,
        "FormatParser< tremolo >::processNeighborInformation() - global id "
        +toString(*iter)+" unknown in additionalAtomData.");
    TremoloAtomInfoContainer &currentInfo = additionalAtomData[*iter];
    currentInfo.torsion = adaptIdDependentDataString(currentInfo.torsion, idglobalizer);
  }
}