source: src/Parser/PdbParser.cpp@ b0a2e3

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

/*
 * PdbParser.cpp
 *
 *  Created on: Aug 17, 2010
 *      Author: heber
 */

// 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 "Descriptors/AtomIdDescriptor.hpp"
#include "World.hpp"
#include "atom.hpp"
#include "bond.hpp"
#include "element.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Parser/PdbParser.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

#include <map>
#include <vector>

#include <iostream>
#include <iomanip>

using namespace std;

/**
 * Constructor.
 */
PdbParser::PdbParser() {
  knownTokens["ATOM"] = PdbKey::Atom;
  knownTokens["HETATM"] = PdbKey::Atom;
  knownTokens["TER"] = PdbKey::Filler;
  knownTokens["END"] = PdbKey::EndOfTimestep;
  knownTokens["CONECT"] = PdbKey::Connect;
  knownTokens["REMARK"] = PdbKey::Remark;
  knownTokens[""] = PdbKey::EndOfTimestep;

  // argh, why can't just PdbKey::X+(size_t)i
  PositionEnumMap[0] = PdbKey::X;
  PositionEnumMap[1] = PdbKey::Y;
  PositionEnumMap[2] = PdbKey::Z;
}

/**
 * Destructor.
 */
PdbParser::~PdbParser() {
  PdbAtomInfoContainer::clearknownDataKeys();
  additionalAtomData.clear();
  atomIdMap.clear();
}


/** Parses the initial word of the given \a line and returns the token type.
 *
 * @param line line to scan
 * @return token type
 */
enum PdbKey::KnownTokens PdbParser::getToken(string &line)
{
  // look for first space
  const size_t space_location = line.find(' ');
  const size_t tab_location = line.find('\t');
  size_t location = space_location < tab_location ? space_location : tab_location;
  string token;
  if (location != string::npos) {
    //DoLog(1) && (Log() << Verbose(1) << "Found space at position " << space_location << std::endl);
    token = line.substr(0,space_location);
  } else {
    token = line;
  }

  //DoLog(1) && (Log() << Verbose(1) << "Token is " << token << std::endl);
  if (knownTokens.count(token) == 0)
    return PdbKey::NoToken;
  else
    return knownTokens[token];

  return PdbKey::NoToken;
}

/**
 * Loads atoms from a PDB-formatted file.
 *
 * \param PDB file
 */
void PdbParser::load(istream* file) {
  string line;
  size_t linecount  = 0;
  enum PdbKey::KnownTokens token;

  // reset atomIdMap for this file (to correctly parse CONECT entries)
  atomIdMap.clear();

  bool NotEndOfFile = true;
  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  unsigned int step = 0;
  // 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 (NotEndOfFile) {
    bool NotEndOfTimestep = true;
    while (NotEndOfTimestep && NotEndOfFile) {
      std::getline(*file, line, '\n');
      if (!line.empty()) {
        // extract first token
        token = getToken(line);
        switch (token) {
          case PdbKey::Atom:
            LOG(3,"INFO: Parsing ATOM entry for time step " << step << ".");
            readAtomDataLine(step, line, newmol);
            break;
          case PdbKey::Remark:
            LOG(3,"INFO: Parsing REM entry for time step " << step << ".");
            break;
          case PdbKey::Connect:
            LOG(3,"INFO: Parsing CONECT entry for time step " << step << ".");
            readNeighbors(step, line);
            break;
          case PdbKey::Filler:
            LOG(3,"INFO: Stumbled upon Filler entry for time step " << step << ".");
            break;
          case PdbKey::EndOfTimestep:
            LOG(3,"INFO: Parsing END entry or empty line for time step " << step << ".");
            NotEndOfTimestep = false;
            break;
          default:
            // TODO: put a throw here
            DoeLog(2) && (eLog() << Verbose(2) << "Unknown token: '" << line << "' for time step " << step << "." << std::endl);
            //ASSERT(0, "PdbParser::load() - Unknown token in line "+toString(linecount)+": "+line+".");
            break;
        }
      }
      NotEndOfFile = NotEndOfFile && (file->good());
      linecount++;
    }
    ++step;
  }
  BOOST_FOREACH(atom *_atom, World::getInstance().getAllAtoms()) {
    for (size_t i=0; i < _atom->getTrajectorySize(); ++i) {
      LOG(2, "INFO: Atom " << _atom->getName() << " is at "
          << _atom->getPositionAtStep(i) << " at time step " << i << ".");
    }
  }
}

/**
 * Saves the \a atoms into as a PDB file.
 *
 * \param file where to save the state
 * \param atoms atoms to store
 */
void PdbParser::save(ostream* file, const std::vector<atom *> &AtomList)
{
  DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl);

  // check for maximum number of time steps
  size_t max_timesteps = 0;
  BOOST_FOREACH(atom *_atom, World::getInstance().getAllAtoms()) {
    for (size_t i=0; i < _atom->getTrajectorySize(); ++i) {
      LOG(2, "INFO: Atom " << _atom->getName() << " is at "
          << _atom->getPositionAtStep(i) << " at time step " << i << ".");
    }
    if (_atom->getTrajectorySize() > max_timesteps)
      max_timesteps = _atom->getTrajectorySize();
  }
  LOG(2,"INFO: Found a maximum of " << max_timesteps << " time steps to store.");

  // re-distribute serials
  // (new atoms might have been added)
  // (serials must be consistent over time steps)
  atomIdMap.clear();
  int AtomNo = 1; // serial number starts at 1 in pdb
  for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
    PdbAtomInfoContainer &atomInfo = getadditionalAtomData(*atomIt);
    setSerial((*atomIt)->getId(), AtomNo);
    atomInfo.set(PdbKey::serial, toString(AtomNo));
    AtomNo++;
  }

  // store all time steps
  for (size_t step = 0; step < max_timesteps; ++step) {
    {
      // add initial remark
      *file << "REMARK created by molecuilder on ";
      time_t now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
      // ctime ends in \n\0, we have to cut away the newline
      std::string time(ctime(&now));
      size_t pos = time.find('\n');
      if (pos != 0)
        *file << time.substr(0,pos);
      else
        *file << time;
      *file << ", time step " << step;
      *file << endl;
    }

    {
      std::map<size_t,size_t> MolIdMap;
      size_t MolNo = 1;  // residue number starts at 1 in pdb
      for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
        const molecule *mol = (*atomIt)->getMolecule();
        if ((mol != NULL) && (MolIdMap.find(mol->getId()) == MolIdMap.end())) {
          MolIdMap[mol->getId()] = MolNo++;
        }
      }
      const size_t MaxMol = MolNo;

      // have a count per element and per molecule (0 is for all homeless atoms)
      std::vector<int> **elementNo = new std::vector<int>*[MaxMol];
      for (size_t i = 0; i < MaxMol; ++i)
        elementNo[i] = new std::vector<int>(MAX_ELEMENTS,1);
      char name[MAXSTRINGSIZE];
      std::string ResidueName;

      // write ATOMs
      for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
        PdbAtomInfoContainer &atomInfo = getadditionalAtomData(*atomIt);
        // gather info about residue
        const molecule *mol = (*atomIt)->getMolecule();
        if (mol == NULL) {
          MolNo = 0;
          atomInfo.set(PdbKey::resSeq, "0");
        } else {
          ASSERT(MolIdMap.find(mol->getId()) != MolIdMap.end(),
              "PdbParser::save() - Mol id "+toString(mol->getId())+" not present despite we set it?!");
          MolNo = MolIdMap[mol->getId()];
          atomInfo.set(PdbKey::resSeq, toString(MolIdMap[mol->getId()]));
          if (atomInfo.get<std::string>(PdbKey::resName) == "-")
            atomInfo.set(PdbKey::resName, mol->getName().substr(0,3));
        }
        // get info about atom
        const size_t  Z = (*atomIt)->getType()->getAtomicNumber();
        if (atomInfo.get<std::string>(PdbKey::name) == "-") {  // if no name set, give it a new name
          sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), (*elementNo[MolNo])[Z]);
          (*elementNo[MolNo])[Z] = ((*elementNo[MolNo])[Z]+1) % 100;   // confine to two digits
          atomInfo.set(PdbKey::name, name);
        }
        // set position
        for (size_t i=0; i<NDIM;++i) {
          stringstream position;
          position << setw(8) << fixed << setprecision(3) << (*atomIt)->getPositionAtStep(step).at(i);
          atomInfo.set(PositionEnumMap[i], position.str());
        }
        // change element and charge if changed
        if (atomInfo.get<std::string>(PdbKey::element) != (*atomIt)->getType()->getSymbol())
          atomInfo.set(PdbKey::element, (*atomIt)->getType()->getSymbol());

        // finally save the line
        saveLine(file, atomInfo);
      }
      for (size_t i = 0; i < MaxMol; ++i)
        delete elementNo[i];
      delete elementNo;

      // write CONECTs
      for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
        writeNeighbors(file, 4, *atomIt);
      }
    }
    // END
    *file << "END" << endl;
  }

}

/** Checks whether there is an entry for the given atom's \a _id.
 *
 * @param _id atom's id we wish to check on
 * @return true - entry present, false - only for atom's father or no entry
 */
bool PdbParser::isPresentadditionalAtomData(unsigned int _id)
{
  return (additionalAtomData.find(_id) != additionalAtomData.end());
}


/** Either returns reference to present entry or creates new with default values.
 *
 * @param _atom atom whose entry we desire
 * @return
 */
PdbAtomInfoContainer& PdbParser::getadditionalAtomData(atom *_atom)
{
  if (additionalAtomData.find(_atom->getId()) != additionalAtomData.end()) {
  } else if (additionalAtomData.find(_atom->father->getId()) != additionalAtomData.end()) {
    // use info from direct father
    additionalAtomData[_atom->getId()] = additionalAtomData[_atom->father->getId()];
  } else if (additionalAtomData.find(_atom->GetTrueFather()->getId()) != additionalAtomData.end()) {
    // use info from topmost father
    additionalAtomData[_atom->getId()] = additionalAtomData[_atom->GetTrueFather()->getId()];
  } else {
    // create new entry use default values if nothing else is known
    additionalAtomData[_atom->getId()] = defaultAdditionalData;
  }
  return additionalAtomData[_atom->getId()];
}

/**
 * Writes one line of PDB-formatted data to the provided stream.
 *
 * \param stream where to write the line to
 * \param *currentAtom the atom of which information should be written
 * \param AtomNo serial number of atom
 * \param *name name of atom, i.e. H01
 * \param ResidueName Name of molecule
 * \param ResidueNo number of residue
 */
void PdbParser::saveLine(
    ostream* file,
    const PdbAtomInfoContainer &atomInfo)
{
  *file << setfill(' ') << left << setw(6)
      << atomInfo.get<std::string>(PdbKey::token);
  *file << setfill(' ') << right << setw(5)
      << atomInfo.get<int>(PdbKey::serial); /* atom serial number */
  *file << " "; /* char 12 is empty */
  *file << setfill(' ') << left << setw(4)
      << atomInfo.get<std::string>(PdbKey::name);  /* atom name */
  *file << setfill(' ') << left << setw(1)
      << atomInfo.get<std::string>(PdbKey::altLoc); /* alternate location/conformation */
  *file << setfill(' ') << left << setw(3)
      << atomInfo.get<std::string>(PdbKey::resName);  /* residue name */
  *file << " "; /* char 21 is empty */
  *file << setfill(' ') << left << setw(1)
      << atomInfo.get<std::string>(PdbKey::chainID); /* chain identifier */
  *file << setfill(' ') << left << setw(4)
      << atomInfo.get<int>(PdbKey::resSeq); /* residue sequence number */
  *file << setfill(' ') << left << setw(1)
      << atomInfo.get<std::string>(PdbKey::iCode); /* iCode */
  *file << "   "; /* char 28-30 are empty */
  // have the following operate on stringstreams such that format specifiers
  // only act on these
  for (size_t i=0;i<NDIM;++i) {
    stringstream position;
    position << fixed << setprecision(3) << showpoint
        << atomInfo.get<double>(PositionEnumMap[i]);
    *file << setfill(' ') << right << setw(8) << position.str();
  }
  {
    stringstream occupancy;
    occupancy << fixed << setprecision(2) << showpoint
        << atomInfo.get<double>(PdbKey::occupancy); /* occupancy */
    *file << setfill(' ') << right << setw(6) << occupancy.str();
  }
  {
    stringstream tempFactor;
    tempFactor << fixed << setprecision(2) << showpoint
        << atomInfo.get<double>(PdbKey::tempFactor); /* temperature factor */
    *file << setfill(' ') << right << setw(6) << tempFactor.str();
  }
  *file << "          "; /* char 68-76 are empty */
  *file << setfill(' ') << right << setw(2) << atomInfo.get<std::string>(PdbKey::element); /* element */
  *file << setfill(' ') << right << setw(2) << atomInfo.get<int>(PdbKey::charge); /* charge */

  *file << endl;
}

/**
 * Writes the neighbor information of one atom to the provided stream.
 *
 * Note that ListOfBonds of WorldTime::CurrentTime is used.
 *
 * \param *file  where to write neighbor information to
 * \param MaxnumberOfNeighbors of neighbors
 * \param *currentAtom to the atom of which to take the neighbor information
 */
void PdbParser::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) {
  int MaxNo = MaxnumberOfNeighbors;
  const BondList & ListOfBonds = currentAtom->getListOfBonds();
  if (!ListOfBonds.empty()) {
    for(BondList::const_iterator currentBond = ListOfBonds.begin(); currentBond != ListOfBonds.end(); ++currentBond) {
      if (MaxNo >= MaxnumberOfNeighbors) {
        *file << "CONECT";
        *file << setw(5) << getSerial(currentAtom->getId());
        MaxNo = 0;
      }
      *file << setw(5) << getSerial((*currentBond)->GetOtherAtom(currentAtom)->getId());
      MaxNo++;
      if (MaxNo == MaxnumberOfNeighbors)
        *file << "\n";
    }
    if (MaxNo != MaxnumberOfNeighbors)
      *file << "\n";
  }
}

/** Retrieves a value from PdbParser::atomIdMap.
 * \param atomid key
 * \return value
 */
size_t PdbParser::getSerial(const size_t atomid) const
{
  ASSERT(atomIdMap.find(atomid) != atomIdMap.end(), "PdbParser::getAtomId: atomid not present in Map.");
  return (atomIdMap.find(atomid)->second);
}

/** Sets an entry in PdbParser::atomIdMap.
 * \param localatomid key
 * \param atomid value
 * \return true - key not present, false - value present
 */
void PdbParser::setSerial(const size_t localatomid, const size_t atomid)
{
  pair<std::map<size_t,size_t>::iterator, bool > inserter;
//  DoLog(1) && (Log() << Verbose(1) << "PdbParser::setAtomId() - Inserting ("
//      << localatomid << " -> " << atomid << ")." << std::endl);
  inserter = atomIdMap.insert( make_pair(localatomid, atomid) );
  ASSERT(inserter.second, "PdbParser::setAtomId: atomId already present in Map.");
}

/** Either returns present atom with given id or a newly created one.
 *
 * @param id_string
 * @return
 */
atom* PdbParser::getAtomToParse(std::string id_string) const
{
  // get the local ID
  ConvertTo<int> toInt;
  unsigned int AtomID = toInt(id_string);
  LOG(4, "INFO: Local id is "+toString(AtomID)+".");
  // get the atomic ID if present
  atom* newAtom = NULL;
  if (atomIdMap.count((size_t)AtomID)) {
    std::map<size_t, size_t>::const_iterator iter = atomIdMap.find(AtomID);
    AtomID = iter->second;
    LOG(4, "INFO: Global id present as " << AtomID << ".");
    // check if atom exists
    newAtom = World::getInstance().getAtom(AtomById(AtomID));
    LOG(5, "INFO: Listing all present atoms with id.");
    BOOST_FOREACH(atom *_atom, World::getInstance().getAllAtoms())
      LOG(5, "INFO: " << *_atom << " with id " << _atom->getId());
  }
  // if not exists, create
  if (newAtom == NULL) {
    newAtom = World::getInstance().createAtom();
    LOG(4, "INFO: No association to global id present, creating atom.");
  } else {
    LOG(4, "INFO: Existing atom found: " << *newAtom << ".");
  }
  return newAtom;
}

void PdbParser::readPdbAtomInfoContainer(PdbAtomInfoContainer &atomInfo, std::string &line) const
{
  LOG(4,"INFO: Parsing token from "+line.substr(0,6)+".");
  atomInfo.set(PdbKey::token, line.substr(0,6));
  LOG(4,"INFO: Parsing serial from "+line.substr(6,5)+".");
  atomInfo.set(PdbKey::serial, line.substr(6,5));

  LOG(4,"INFO: Parsing name from "+line.substr(12,4)+".");
  atomInfo.set(PdbKey::name, line.substr(12,4));
  LOG(4,"INFO: Parsing altLoc from "+line.substr(16,1)+".");
  atomInfo.set(PdbKey::altLoc, line.substr(16,1));
  LOG(4,"INFO: Parsing resName from "+line.substr(17,3)+".");
  atomInfo.set(PdbKey::resName, line.substr(17,3));
  LOG(4,"INFO: Parsing chainID from "+line.substr(21,1)+".");
  atomInfo.set(PdbKey::chainID, line.substr(21,1));
  LOG(4,"INFO: Parsing resSeq from "+line.substr(22,4)+".");
  atomInfo.set(PdbKey::resSeq, line.substr(22,4));
  LOG(4,"INFO: Parsing iCode from "+line.substr(26,1)+".");
  atomInfo.set(PdbKey::iCode, line.substr(26,1));

  LOG(4,"INFO: Parsing occupancy from "+line.substr(54,6)+".");
  atomInfo.set(PdbKey::occupancy, line.substr(54,6));
  LOG(4,"INFO: Parsing tempFactor from "+line.substr(60,6)+".");
  atomInfo.set(PdbKey::tempFactor, line.substr(60,6));
  LOG(4,"INFO: Parsing charge from "+line.substr(78,2)+".");
  atomInfo.set(PdbKey::charge, line.substr(78,2));
  LOG(4,"INFO: Parsing element from "+line.substr(76,2)+".");
  atomInfo.set(PdbKey::element, line.substr(76,2));
}

/** Parse an ATOM line from a PDB file.
 *
 * Reads one data line of a pdstatus file and interprets it according to the
 * specifications of the PDB 3.2 format: http://www.wwpdb.org/docs.html
 *
 *  A new atom is created and filled with available information, non-
 *  standard information is placed in additionalAtomData at the atom's id.
 *
 * \param _step time step to use
 * \param line to parse as an atom
 * \param newmol molecule to add parsed atoms to
 */
void PdbParser::readAtomDataLine(const unsigned int _step, std::string &line, molecule *newmol = NULL) {
  vector<string>::iterator it;

  atom* newAtom = getAtomToParse(line.substr(6,5));
  LOG(3,"INFO: Parsing END entry or empty line.");
  bool FirstTimestep = isPresentadditionalAtomData(newAtom->getId()) ? false : true;
  ASSERT((FirstTimestep && (_step == 0)) || (!FirstTimestep && (_step !=0)),
      "PdbParser::readAtomDataLine() - logig mismatch between FirstTimestep and step == 0.");
  if (FirstTimestep) {
    LOG(3,"INFO: Parsing new atom.");
  } else {
    LOG(3,"INFO: Parsing present atom "+toString(*newAtom)+".");
  }
  PdbAtomInfoContainer &atomInfo = getadditionalAtomData(newAtom);
  LOG(4,"INFO: Information in container is "+toString(atomInfo)+".");

  string word;
  ConvertTo<size_t> toSize_t;

  // assign highest+1 instead, but then beware of CONECT entries! Another map needed!
//  if (!Inserter.second) {
//    const size_t id = (*SerialSet.rbegin())+1;
//    SerialSet.insert(id);
//    atomInfo.set(PdbKey::serial, toString(id));
//    DoeLog(2) && (eLog() << Verbose(2)
//        << "Serial " << atomInfo.get<std::string>(PdbKey::serial) << " already present, "
//        << "assigning " << toString(id) << " instead." << std::endl);
//  }

  // check whether serial exists, if so, assign next available

//  DoLog(2) && (Log() << Verbose(2) << "Split line:"
//      << line.substr(6,5) << "|"
//      << line.substr(12,4) << "|"
//      << line.substr(16,1) << "|"
//      << line.substr(17,3) << "|"
//      << line.substr(21,1) << "|"
//      << line.substr(22,4) << "|"
//      << line.substr(26,1) << "|"
//      << line.substr(30,8) << "|"
//      << line.substr(38,8) << "|"
//      << line.substr(46,8) << "|"
//      << line.substr(54,6) << "|"
//      << line.substr(60,6) << "|"
//      << line.substr(76,2) << "|"
//      << line.substr(78,2) << std::endl);

  if (FirstTimestep) {
    // first time step
    // then fill info container
    readPdbAtomInfoContainer(atomInfo, line);
    // set the serial
    std::pair< std::set<size_t>::const_iterator, bool> Inserter =
      SerialSet.insert(toSize_t(atomInfo.get<std::string>(PdbKey::serial)));
    ASSERT(Inserter.second,
        "PdbParser::readAtomDataLine() - ATOM contains entry with serial "
        +atomInfo.get<std::string>(PdbKey::serial)+" already present!");
    setSerial(toSize_t(atomInfo.get<std::string>(PdbKey::serial)), newAtom->getId());
    // set position
    Vector tempVector;
    LOG(4,"INFO: Parsing position from ("
        +line.substr(30,8)+","
        +line.substr(38,8)+","
        +line.substr(46,8)+").");
    PdbAtomInfoContainer::ScanKey(tempVector[0], line.substr(30,8));
    PdbAtomInfoContainer::ScanKey(tempVector[1], line.substr(38,8));
    PdbAtomInfoContainer::ScanKey(tempVector[2], line.substr(46,8));
    newAtom->setPosition(tempVector);
    // set element
    const element *elem = World::getInstance().getPeriode()
        ->FindElement(atomInfo.get<std::string>(PdbKey::element));
    ASSERT(elem != NULL,
        "PdbParser::readAtomDataLine() - element "+atomInfo.get<std::string>(PdbKey::element)+" is unknown!");
    newAtom->setType(elem);

    if (newmol != NULL)
      newmol->AddAtom(newAtom);
  } else {
    // not first time step
    // then parse into different container
    PdbAtomInfoContainer consistencyInfo;
    readPdbAtomInfoContainer(consistencyInfo, line);
    // then check additional info for consistency
    ASSERT(atomInfo.get<std::string>(PdbKey::token) == consistencyInfo.get<std::string>(PdbKey::token),
        "PdbParser::readAtomDataLine() - difference in token on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::name) == consistencyInfo.get<std::string>(PdbKey::name),
        "PdbParser::readAtomDataLine() - difference in name on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+":"
        +atomInfo.get<std::string>(PdbKey::name)+"!="
        +consistencyInfo.get<std::string>(PdbKey::name)+".");
    ASSERT(atomInfo.get<std::string>(PdbKey::altLoc) == consistencyInfo.get<std::string>(PdbKey::altLoc),
        "PdbParser::readAtomDataLine() - difference in altLoc on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::resName) == consistencyInfo.get<std::string>(PdbKey::resName),
        "PdbParser::readAtomDataLine() - difference in resName on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::chainID) == consistencyInfo.get<std::string>(PdbKey::chainID),
        "PdbParser::readAtomDataLine() - difference in chainID on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::resSeq) == consistencyInfo.get<std::string>(PdbKey::resSeq),
        "PdbParser::readAtomDataLine() - difference in resSeq on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::iCode) == consistencyInfo.get<std::string>(PdbKey::iCode),
        "PdbParser::readAtomDataLine() - difference in iCode on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::occupancy) == consistencyInfo.get<std::string>(PdbKey::occupancy),
        "PdbParser::readAtomDataLine() - difference in occupancy on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::tempFactor) == consistencyInfo.get<std::string>(PdbKey::tempFactor),
        "PdbParser::readAtomDataLine() - difference in tempFactor on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::charge) == consistencyInfo.get<std::string>(PdbKey::charge),
        "PdbParser::readAtomDataLine() - difference in charge on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    ASSERT(atomInfo.get<std::string>(PdbKey::element) == consistencyInfo.get<std::string>(PdbKey::element),
        "PdbParser::readAtomDataLine() - difference in element on multiple time step for atom with id "
        +atomInfo.get<std::string>(PdbKey::serial)+"!");
    // and parse in trajectory
    Vector tempVector;
    LOG(4,"INFO: Parsing trajectory position from ("
        +line.substr(30,8)+","
        +line.substr(38,8)+","
        +line.substr(46,8)+").");
    PdbAtomInfoContainer::ScanKey(tempVector[0], line.substr(30,8));
    PdbAtomInfoContainer::ScanKey(tempVector[1], line.substr(38,8));
    PdbAtomInfoContainer::ScanKey(tempVector[2], line.substr(46,8));
    LOG(4,"INFO: Adding trajectory point to time step "+toString(_step)+".");
    // and set position at new time step
    newAtom->setPositionAtStep(_step, tempVector);
  }


//  printAtomInfo(newAtom);
}

/** Prints all PDB-specific information known about an atom.
 *
 */
void PdbParser::printAtomInfo(const atom * const newAtom) const
{
  const PdbAtomInfoContainer &atomInfo = additionalAtomData.at(newAtom->getId()); // operator[] const does not exist

  DoLog(1) && (Log() << Verbose(1) << "We know about atom " << newAtom->getId() << ":" << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\ttoken is " << atomInfo.get<std::string>(PdbKey::token) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tserial is " << atomInfo.get<int>(PdbKey::serial) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tname is " << atomInfo.get<std::string>(PdbKey::name) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\taltLoc is " << atomInfo.get<std::string>(PdbKey::altLoc) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tresName is " << atomInfo.get<std::string>(PdbKey::resName) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tchainID is " << atomInfo.get<std::string>(PdbKey::chainID) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tresSeq is " << atomInfo.get<int>(PdbKey::resSeq) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tiCode is " << atomInfo.get<std::string>(PdbKey::iCode) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tX is " << atomInfo.get<double>(PdbKey::X) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tY is " << atomInfo.get<double>(PdbKey::Y) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tZ is " << atomInfo.get<double>(PdbKey::Z) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\toccupancy is " << atomInfo.get<double>(PdbKey::occupancy) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\ttempFactor is " << atomInfo.get<double>(PdbKey::tempFactor) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\telement is '" << *(newAtom->getType()) << "'" << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tcharge is " << atomInfo.get<int>(PdbKey::charge) << std::endl);
}

/**
 * Reads neighbor information for one atom from the input.
 *
 * \param _step time step to use
 * \param line to parse as an atom
 */
void PdbParser::readNeighbors(const unsigned int _step, std::string &line)
{
  const size_t length = line.length();
  std::list<size_t> ListOfNeighbors;
  ConvertTo<size_t> toSize_t;

  // obtain neighbours
  // show split line for debugging
  string output;
  ASSERT(length >=16,
      "PdbParser::readNeighbors() - CONECT entry has not enough entries: "+line+"!");
//  output = "Split line:|";
//  output += line.substr(6,5) + "|";
  const size_t id = toSize_t(line.substr(6,5));
  for (size_t index = 11; index <= 26; index+=5) {
    if (index+5 <= length) {
//      output += line.substr(index,5) + "|";
      const size_t otherid = toSize_t(line.substr(index,5));
      ListOfNeighbors.push_back(otherid);
    } else  {
      break;
    }
  }
//  DoLog(2) && (Log() << Verbose(2) << output << std::endl);

  // add neighbours
  atom *_atom = World::getInstance().getAtom(AtomById(getSerial(id)));
  for (std::list<size_t>::const_iterator iter = ListOfNeighbors.begin();
      iter != ListOfNeighbors.end();
      ++iter) {
//    DoLog(1) && (Log() << Verbose(1) << "Adding Bond (" << getAtomId(id) << "," << getAtomId(*iter) << ")" << std::endl);
    atom * const _Otheratom = World::getInstance().getAtom(AtomById(getSerial(*iter)));
    _atom->addBond(_step, _Otheratom);
  }
}

/**
 * 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
 */
//string PdbParser::adaptIdDependentDataString(string data) {
//  // there might be no IDs
//  if (data == "-") {
//    return "-";
//  }
//
//  char separator;
//  int id;
//  stringstream line, result;
//  line << data;
//
//  line >> id;
//  result << atomIdMap[id];
//  while (line.good()) {
//    line >> separator >> id;
//    result << separator << atomIdMap[id];
//  }
//
//  return result.str();
//  return "";
//}


bool PdbParser::operator==(const PdbParser& b) const
{
  bool status = true;
  World::AtomComposite atoms = World::getInstance().getAllAtoms();
  for (World::AtomComposite::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    if ((additionalAtomData.find((*iter)->getId()) != additionalAtomData.end())
        && (b.additionalAtomData.find((*iter)->getId()) != b.additionalAtomData.end())) {
      const PdbAtomInfoContainer &atomInfo = additionalAtomData.at((*iter)->getId());
      const PdbAtomInfoContainer &OtheratomInfo = b.additionalAtomData.at((*iter)->getId());

      status = status && (atomInfo.get<std::string>(PdbKey::serial) == OtheratomInfo.get<std::string>(PdbKey::serial));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in serials!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::name) == OtheratomInfo.get<std::string>(PdbKey::name));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in names!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::altLoc) == OtheratomInfo.get<std::string>(PdbKey::altLoc));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in altLocs!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::resName) == OtheratomInfo.get<std::string>(PdbKey::resName));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resNames!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::chainID) == OtheratomInfo.get<std::string>(PdbKey::chainID));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in chainIDs!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::resSeq) == OtheratomInfo.get<std::string>(PdbKey::resSeq));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resSeqs!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::iCode) == OtheratomInfo.get<std::string>(PdbKey::iCode));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in iCodes!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::occupancy) == OtheratomInfo.get<std::string>(PdbKey::occupancy));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in occupancies!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::tempFactor) == OtheratomInfo.get<std::string>(PdbKey::tempFactor));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in tempFactors!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::charge) == OtheratomInfo.get<std::string>(PdbKey::charge));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in charges!" << std::endl);
    }
  }

  return status;
}