source: src/Parser/TremoloParser.cpp@ c3ecde

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * 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 "LinearAlgebra/RealSpaceMatrix.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "World.hpp"
#include "WorldTime.hpp"


#include <algorithm>
#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)
{
  createKnownTypesByIdentity();

  // 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;
  // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
  World::getInstance().getMolecules()->insert(newmol);
  while (file->good()) {
    std::getline(*file, line, '\n');
    // we only parse in the first ATOMDATA line
    if (usedFields_load.empty()) {
      location = line.find("ATOMDATA", 0);
      if (location != string::npos) {
       parseAtomDataKeysLine(line, location + 8, usedFields_load);
      }
    }
    if (line.length() > 0 && line.at(0) != '#') {
      readAtomDataLine(line, newmol);
    }
  }
  LOG(3, "DEBUG: Local usedFields is: " << usedFields_load);

  // 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();
  adaptTorsion();

  // 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<atom *> &AtomList) {
  LOG(0, "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, "DEBUG: Global (with unique entries) usedFields_save is: " << usedFields_save);

  // distribute ids continuously
  distributeContinuousIds(AtomList);

  // store atomdata
  save_AtomDataLine(file);

  // store box
  save_BoxLine(file);

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

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<atom *> &AtomList)
{
  resetIdAssociations();
  atomId_t lastid = 0;
  for (std::vector<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);
}

/**
 * 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* currentAtom)
{
  TremoloKey::atomDataKey currentField;

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

  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->getPosition());
        *file << currentAtom->at(0) << "\t";
        *file << currentAtom->at(1) << "\t";
        *file << currentAtom->at(2) << "\t";
        break;
      case TremoloKey::u :
        // for the moment, assume there are always three dimensions
        LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicVelocity());
        *file << currentAtom->getAtomicVelocity()[0] << "\t";
        *file << currentAtom->getAtomicVelocity()[1] << "\t";
        *file << currentAtom->getAtomicVelocity()[2] << "\t";
        break;
      case TremoloKey::F :
        // for the moment, assume there are always three dimensions
        LOG(3, "Writing for type " << knownKeyNames[currentField] << ": " << currentAtom->getAtomicForce());
        *file << currentAtom->getAtomicForce()[0] << "\t";
        *file << currentAtom->getAtomicForce()[1] << "\t";
        *file << currentAtom->getAtomicForce()[2] << "\t";
        break;
      case TremoloKey::type :
        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) << "\t";
          } else {
            LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value: " << currentAtom->getType()->getSymbol());
            *file << currentAtom->getType()->getSymbol() << "\t";
          }
        } 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) << "\t";
          } else {
            LOG(3, "Writing for type " << knownKeyNames[currentField] << " default value from father: " << currentAtom->GetTrueFather()->getType()->getSymbol());
            *file << currentAtom->GetTrueFather()->getType()->getSymbol() << "\t";
          }
        } else {
          LOG(3, "Writing for type " << knownKeyNames[currentField] << " its default value: " << currentAtom->getType()->getSymbol());
          *file << currentAtom->getType()->getSymbol() << "\t";
        }
        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::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()+1;
        } 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.) {
          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] << " AtomInfo::charge : " << currentAtom->getCharge());
            *file << currentAtom->getCharge();
          }
        } else {
          LOG(3, "Writing for type " << knownKeyNames[currentField] << " AtomInfo::charge : " << currentAtom->getCharge());
          *file << currentAtom->getCharge();
        }
        *file << "\t";
        break;
      default :
        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] << " the default: " << defaultAdditionalData.get(currentField));
          *file << defaultAdditionalData.get(currentField);
        }
        *file << "\t";
        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* 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);
  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);
  }
  //LOG(1, "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
 */
void FormatParser< tremolo >::readAtomDataLine(const std::string &line, molecule *newmol = NULL) {
  std::stringstream lineStream;
  atom* newAtom = World::getInstance().createAtom();
  const atomId_t atomid = newAtom->getId();
  additionalAtomData[atomid] = TremoloAtomInfoContainer(); // fill with default values
  TremoloAtomInfoContainer *atomInfo = &additionalAtomData[atomid];
  TremoloKey::atomDataKey currentField;
  ConvertTo<double> toDouble;
  ConvertTo<int> toInt;
  Vector tempVector;

  // 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::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];
    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 " << *tok_iter);
          newAtom->set(i, toDouble(*tok_iter));
          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 " << *tok_iter);
          tempVector[i] = toDouble(*tok_iter);
          tok_iter++;
        }
        newAtom->setAtomicVelocity(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 " << *tok_iter);
          tempVector[i] = toDouble(*tok_iter);
          tok_iter++;
        }
        newAtom->setAtomicForce(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 " << *tok_iter);
        std::string element;
        try {
          element = knownTypes.getType(*tok_iter);
        } catch(IllegalParserKeyException) {
          // clean up
          World::getInstance().destroyAtom(newAtom);
          // give an error
          ELOG(0, "TremoloParser: I do not understand the element token " << *tok_iter << ".");
        }
        // put type name into container for later use
        atomInfo->set(currentField, *tok_iter);
        LOG(4, "INFO: Parsing element " << (*tok_iter) << " as " << element << " according to KnownTypes.");
        tok_iter++;
        newAtom->setType(World::getInstance().getPeriode()->FindElement(element));
        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 " << *tok_iter);
        associateLocaltoGlobalId(toInt(*tok_iter), atomid);
        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+"!");
          LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
          lineStream << *tok_iter << "\t";
          tok_iter++;
        }
        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+"!");
        LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
        atomInfo->set(currentField, *tok_iter);
        newAtom->setCharge(boost::lexical_cast<double>(*tok_iter));
        tok_iter++;
        break;
      default :
        ASSERT(tok_iter != tokens.end(), "FormatParser< tremolo >::readAtomDataLine() - no value for "+keyName+"!");
        LOG(4, "INFO: Parsing key " << keyName << " with next token " << *tok_iter);
        atomInfo->set(currentField, *tok_iter);
        tok_iter++;
        break;
    }
  }
  LOG(3, "INFO: Parsed atom " << atomid << ".");
  if (newmol != NULL)
    newmol->AddAtom(newAtom);
}

bool FormatParser< tremolo >::saveAtomsInExttypes(std::ostream &output, const std::vector<atom*> &atoms, const int id) const
{
  bool status = true;
  // parse the file
  for (std::vector<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
 *
 * \return input string with modified atom IDs
 */
std::string FormatParser< tremolo >::adaptIdDependentDataString(std::string data) {
  // there might be no IDs
  if (data == "-") {
    return "-";
  }

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

  line >> id;
  result << getGlobalId(id);
  while (line.good()) {
    line >> separator >> id;
    result << separator << getGlobalId(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.
 */
void FormatParser< tremolo >::adaptImprData() {
  if (!isUsedField(usedFields_load, "imprData")) {
    return;
  }

  for(std::map<const atomId_t, TremoloAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
    currentInfo != additionalAtomData.end(); currentInfo++
  ) {
    currentInfo->second.imprData = adaptIdDependentDataString(currentInfo->second.imprData);
  }
}

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

  for(std::map<const atomId_t, TremoloAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
    currentInfo != additionalAtomData.end(); currentInfo++
  ) {
    currentInfo->second.torsion = adaptIdDependentDataString(currentInfo->second.torsion);
  }
}