source: src/Parser/PdbParser.cpp@ 286a967

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. 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/>.
 */

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

#include <map>
#include <vector>

#include <iostream>
#include <iomanip>

using namespace std;

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

/**
 * Constructor.
 */
FormatParser< pdb >::FormatParser() :
  FormatParser_common(NULL)
{
  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.
 */
FormatParser< pdb >::~FormatParser() 
{
  PdbAtomInfoContainer::clearknownDataKeys();
  additionalAtomData.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 FormatParser< pdb >::getToken(std::string &line)
{
  // look for first space
  std::string token = line.substr(0,6);
  const size_t space_location = token.find(' ');
  const size_t tab_location = token.find('\t');
  size_t location = space_location < tab_location ? space_location : tab_location;
  if (location != string::npos) {
    //LOG(1, "Found space at position " << space_location);
    token = token.substr(0,space_location);
  }

  //LOG(1, "Token is " << token);
  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 FormatParser< pdb >::load(istream* file) {
  string line;
  size_t linecount  = 0;
  enum PdbKey::KnownTokens token;

  // reset id maps for this file (to correctly parse CONECT entries)
  resetIdAssociations();

  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(1,"INFO: Parsing END entry or empty line for time step " << step << ".");
            NotEndOfTimestep = false;
            break;
          default:
            // TODO: put a throw here
            ELOG(2, "Unknown token: '" << line << "' for time step " << step << ".");
            //ASSERT(0, "FormatParser< pdb >::load() - Unknown token in line "+toString(linecount)+": "+line+".");
            break;
        }
      }
      NotEndOfFile = NotEndOfFile && (file->good());
      linecount++;
    }
    ++step;
  }
  LOG(4, "INFO: Listing all newly parsed atoms.");
  BOOST_FOREACH(atom *_atom, *newmol)
    LOG(4, "INFO: Atom " << _atom->getName() << " " << *dynamic_cast<AtomInfo *>(_atom) << ".");

  // refresh atom::nr and atom::name
  newmol->getAtomCount();
}

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

  // check for maximum number of time steps
  size_t max_timesteps = 0;
  BOOST_FOREACH(atom *_atom, World::getInstance().getAllAtoms()) {
    LOG(4, "INFO: Atom " << _atom->getName() << " " << *dynamic_cast<AtomInfo *>(_atom) << ".");
    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
  resetIdAssociations();
  // (new atoms might have been added)
  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);
    associateLocaltoGlobalId(AtomNo, (*atomIt)->getId());
    atomInfo.set(PdbKey::serial, toString(AtomNo));
    ++AtomNo;
  }

  // store all time steps (always do first step)
  for (size_t step = 0; (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(),
              "FormatParser< pdb >::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()) {
          std::string symbol = (*atomIt)->getType()->getSymbol();
          if ((symbol[1] >= 'a') && (symbol[1] <= 'z'))
            symbol[1] = (symbol[1] - 'a') + 'A';
          atomInfo.set(PdbKey::element, symbol);
        }

        // 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;
  }

}

/** Add default info, when new atom is added to World.
 *
 * @param id of atom
 */
void FormatParser< pdb >::AtomInserted(atomId_t id)
{
  //LOG(3, "FormatParser< pdb >::AtomInserted() - notified of atom " << id << "'s insertion.");
  ASSERT(!isPresentadditionalAtomData(id),
      "FormatParser< pdb >::AtomInserted() - additionalAtomData already present for newly added atom "
      +toString(id)+".");
  // don't insert here as this is our check whether we are in the first time step
  //additionalAtomData.insert( std::make_pair(id, defaultAdditionalData) );
}

/** Remove additional AtomData info, when atom has been removed from World.
 *
 * @param id of atom
 */
void FormatParser< pdb >::AtomRemoved(atomId_t id)
{
  //LOG(3, "FormatParser< pdb >::AtomRemoved() - notified of atom " << id << "'s removal.");
  std::map<const atomId_t, PdbAtomInfoContainer>::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< pdb >::AtomRemoved() - additionalAtomData is not present for atom "
//      +toString(id)+" to remove.");
  if (iter != additionalAtomData.end()) {
    additionalAtomData.erase(iter);
  }
}


/** 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 FormatParser< pdb >::isPresentadditionalAtomData(const atomId_t _id) const
{
  std::map<const atomId_t, PdbAtomInfoContainer>::const_iterator iter = additionalAtomData.find(_id);
  return (iter != additionalAtomData.end());
}


/** Either returns reference to present entry or creates new with default values.
 *
 * @param _atom atom whose entry we desire
 * @return
 */
PdbAtomInfoContainer& FormatParser< pdb >::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 FormatParser< pdb >::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.
 *
 * Also, we fill up the CONECT line to extend over 80 chars.
 *
 * \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 FormatParser< pdb >::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) {
  int MaxNo = MaxnumberOfNeighbors;
  int charsleft = 80;
  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) << getLocalId(currentAtom->getId());
        charsleft = 80-6-5;
        MaxNo = 0;
      }
      *file << setw(5) << getLocalId((*currentBond)->GetOtherAtom(currentAtom)->getId());
      charsleft -= 5;
      MaxNo++;
      if (MaxNo == MaxnumberOfNeighbors) {
        for (;charsleft > 0; charsleft--)
          *file << ' ';
        *file << "\n";
      }
    }
    if (MaxNo != MaxnumberOfNeighbors) {
      for (;charsleft > 0; charsleft--)
        *file << ' ';
      *file << "\n";
    }
  }
}

/** Either returns present atom with given id or a newly created one.
 *
 * @param id_string
 * @return
 */
atom* FormatParser< pdb >::getAtomToParse(std::string id_string)
{
  // get the local ID
  ConvertTo<int> toInt;
  const unsigned int AtomID_local = toInt(id_string);
  LOG(4, "INFO: Local id is "+toString(AtomID_local)+".");
  // get the atomic ID if present
  atom* newAtom = NULL;
  if (getGlobalId(AtomID_local) != -1) {
    const unsigned int AtomID_global = getGlobalId(AtomID_local);
    LOG(4, "INFO: Global id present as " << AtomID_global << ".");
    // check if atom exists
    newAtom = World::getInstance().getAtom(AtomById(AtomID_global));
    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();
    //const unsigned int AtomID_global = newAtom->getId();
    LOG(4, "INFO: No association to global id present, creating atom.");
  } else {
    LOG(4, "INFO: Existing atom found: " << *newAtom << ".");
  }
  return newAtom;
}

/** read a line starting with key ATOM.
 *
 * We check for line's length and parse only up to this value.
 *
 * @param atomInfo container to put information in
 * @param line line containing key ATOM
 */
void FormatParser< pdb >::readPdbAtomInfoContainer(PdbAtomInfoContainer &atomInfo, std::string &line) const
{
  const size_t length = line.length();
  if (length < 80)
    ELOG(2, "FormatParser< pdb >::readPdbAtomInfoContainer() - pdb is mal-formed, containing less than 80 chars!");
  if (length >= 6) {
    LOG(4,"INFO: Parsing token from "+line.substr(0,6)+".");
    atomInfo.set(PdbKey::token, line.substr(0,6));
  }
  if (length >= 11) {
    LOG(4,"INFO: Parsing serial from "+line.substr(6,5)+".");
    atomInfo.set(PdbKey::serial, line.substr(6,5));
    ASSERT(atomInfo.get<int>(PdbKey::serial) != 0,
        "FormatParser< pdb >::readPdbAtomInfoContainer() - serial 0 is invalid (filler id for conect entries).");
  }

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

  if (length >= 60) {
    LOG(4,"INFO: Parsing occupancy from "+line.substr(54,6)+".");
    atomInfo.set(PdbKey::occupancy, line.substr(54,6));
  }
  if (length >= 66) {
    LOG(4,"INFO: Parsing tempFactor from "+line.substr(60,6)+".");
    atomInfo.set(PdbKey::tempFactor, line.substr(60,6));
  }
  if (length >= 80) {
    LOG(4,"INFO: Parsing charge from "+line.substr(78,2)+".");
    atomInfo.set(PdbKey::charge, line.substr(78,2));
  }
  if (length >= 78) {
    LOG(4,"INFO: Parsing element from "+line.substr(76,2)+".");
    atomInfo.set(PdbKey::element, line.substr(76,2));
  } else {
    LOG(4,"INFO: Trying to parse alternative element from name "+line.substr(12,4)+".");
    atomInfo.set(PdbKey::element, line.substr(12,4));
  }
}

/** 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 FormatParser< pdb >::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)),
      "FormatParser< pdb >::readAtomDataLine() - additionalAtomData present though atom is newly parsed.");
  if (FirstTimestep) {
    LOG(3,"INFO: Parsing new atom "+toString(*newAtom)+" "+toString(newAtom->getId())+".");
  } 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;

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

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

  if (FirstTimestep) {
    // first time step
    // then fill info container
    readPdbAtomInfoContainer(atomInfo, line);
    // associate local with global id
    associateLocaltoGlobalId(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
    std::string value = atomInfo.get<std::string>(PdbKey::element);
    // make second character lower case if not
    if ((value[1] >= 'A') && (value[1] <= 'Z'))
      value[1] = (value[1] - 'A') + 'a';
    const element *elem = World::getInstance().getPeriode()
        ->FindElement(value);
    ASSERT(elem != NULL,
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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),
        "FormatParser< pdb >::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 FormatParser< pdb >::printAtomInfo(const atom * const newAtom) const
{
  const PdbAtomInfoContainer &atomInfo = additionalAtomData.at(newAtom->getId()); // operator[] const does not exist

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

/**
 * Reads neighbor information for one atom from the input.
 *
 * \param _step time step to use
 * \param line to parse as an atom
 */
void FormatParser< pdb >::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,
      "FormatParser< pdb >::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) + "|";
      // search for digits
      int otherid = -1;
      PdbAtomInfoContainer::ScanKey(otherid, line.substr(index,5));
      if (otherid > 0)
        ListOfNeighbors.push_back(otherid);
      else
        ELOG(3, "FormatParser< pdb >::readNeighbors() - discarding CONECT entry with id 0.");
    } else  {
      break;
    }
  }
  LOG(4, output);

  // add neighbours
  atom *_atom = World::getInstance().getAtom(AtomById(getGlobalId(id)));
  LOG(2, "STATUS: Atom " << _atom->getId() << " gets " << ListOfNeighbors.size() << " more neighbours.");
  for (std::list<size_t>::const_iterator iter = ListOfNeighbors.begin();
      iter != ListOfNeighbors.end();
      ++iter) {
    atom * const _Otheratom = World::getInstance().getAtom(AtomById(getGlobalId(*iter)));
    LOG(3, "INFO: Adding Bond (" << *_atom << "," << *_Otheratom << ")");
    _atom->addBond(_step, _Otheratom);
  }
}

bool FormatParser< pdb >::operator==(const FormatParser< pdb >& 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) ELOG(1, "Mismatch in serials!");
      status = status && (atomInfo.get<std::string>(PdbKey::name) == OtheratomInfo.get<std::string>(PdbKey::name));
      if (!status) ELOG(1, "Mismatch in names!");
      status = status && (atomInfo.get<std::string>(PdbKey::altLoc) == OtheratomInfo.get<std::string>(PdbKey::altLoc));
      if (!status) ELOG(1, "Mismatch in altLocs!");
      status = status && (atomInfo.get<std::string>(PdbKey::resName) == OtheratomInfo.get<std::string>(PdbKey::resName));
      if (!status) ELOG(1, "Mismatch in resNames!");
      status = status && (atomInfo.get<std::string>(PdbKey::chainID) == OtheratomInfo.get<std::string>(PdbKey::chainID));
      if (!status) ELOG(1, "Mismatch in chainIDs!");
      status = status && (atomInfo.get<std::string>(PdbKey::resSeq) == OtheratomInfo.get<std::string>(PdbKey::resSeq));
      if (!status) ELOG(1, "Mismatch in resSeqs!");
      status = status && (atomInfo.get<std::string>(PdbKey::iCode) == OtheratomInfo.get<std::string>(PdbKey::iCode));
      if (!status) ELOG(1, "Mismatch in iCodes!");
      status = status && (atomInfo.get<std::string>(PdbKey::occupancy) == OtheratomInfo.get<std::string>(PdbKey::occupancy));
      if (!status) ELOG(1, "Mismatch in occupancies!");
      status = status && (atomInfo.get<std::string>(PdbKey::tempFactor) == OtheratomInfo.get<std::string>(PdbKey::tempFactor));
      if (!status) ELOG(1, "Mismatch in tempFactors!");
      status = status && (atomInfo.get<std::string>(PdbKey::charge) == OtheratomInfo.get<std::string>(PdbKey::charge));
      if (!status) ELOG(1, "Mismatch in charges!");
    }
  }

  return status;
}