source: src/Parser/PdbParser.cpp@ 24f128

/*
 * 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 "Helpers/MemDebug.hpp"

#include "Helpers/Assert.hpp"
#include "Helpers/Log.hpp"
#include "Helpers/toString.hpp"
#include "Helpers/Verbose.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 <map>
#include <vector>

#include <iostream>
#include <iomanip>

using namespace std;

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

/**
 * Destructor.
 */
PdbParser::~PdbParser() {
  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;

  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  bool NotEndOfFile = 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 (NotEndOfFile) {
    std::getline(*file, line, '\n');
    // extract first token
    token = getToken(line);
    //DoLog(1) && (Log() << Verbose(1) << " Recognized token of type : " << token << std::endl);
    switch (token) {
      case PdbKey::Atom:
        readAtomDataLine(line, newmol);
        break;
      case PdbKey::Remark:
        break;
      case PdbKey::Connect:
        readNeighbors(line);
        break;
      case PdbKey::Filler:
        break;
      case PdbKey::EndOfFile:
        NotEndOfFile = false;
        break;
      default:
        // TODO: put a throw here
        DoeLog(2) && (eLog() << Verbose(2) << "Unknown token: '" << line << "'" << std::endl);
        //ASSERT(0, "PdbParser::load() - Unknown token in line "+toString(linecount)+": "+line+".");
        break;
    }
    NotEndOfFile = NotEndOfFile && (file->good());
    linecount++;
  }
}

/**
 * Saves the World's current state into as a PDB file.
 *
 * \param file where to save the state
 */
void PdbParser::save(ostream* file) {
  DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl);

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

  // we distribute new atom numbers, hence clear map beforehand
  atomIdMap.clear();
  {
    vector<atom *> AtomList = World::getInstance().getAllAtoms();

    std::vector<int> elementNo(MAX_ELEMENTS,1);
    char name[MAXSTRINGSIZE];

    // write ATOMs
    int AtomNo = 1; // serial number starts at 1 in pdb
    int MolNo = 1;  // residue number starts at 1 in pdb
    for (vector<atom *>::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
      const size_t  Z = (*atomIt)->getType()->getAtomicNumber();
      sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), elementNo[Z]);
      elementNo[Z] = (elementNo[Z]+1) % 100;   // confine to two digits
      const molecule *mol = (*atomIt)->getMolecule();
      if (mol == NULL) {  // for homeless atoms, MolNo = -1 is reserved
        MolNo = -1;
      } else {
        MolNo = mol->getId();
      }
      saveLine(file, *atomIt, name, AtomNo, MolNo);
      setAtomId((*atomIt)->getId(), AtomNo);
      AtomNo++;
    }

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

  // END
  *file << "END" << endl;
}

/**
 * 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 *name name of atom, i.e. H01
 * \param AtomNo serial number of atom
 * \param ResidueNo number of residue
 */
void PdbParser::saveLine(ostream* file, const atom* currentAtom, const char *name, const int AtomNo, const int ResidueNo) {
  *file << "ATOM ";
  *file << setw(6) << AtomNo; /* atom serial number */
  *file << setw(1) << " ";
  *file << setfill(' ') << left << setw(4) << name << right;  /* atom name */
  *file << setw(1) << " ";
  *file << setfill(' ') << setw(3) << ((currentAtom->getMolecule() != NULL) ? currentAtom->getMolecule()->getName().substr(0,3) : "-");  /* residue name */
  *file << setw(1) << " ";
  *file << setfill(' ') << setw(1) << (char)('a'+(unsigned char)(AtomNo % 26)); /* letter for chain */
  *file << setw(4) << ResidueNo; /* residue sequence number */
  *file << setw(4) << "    ";
  for (int i=0;i<NDIM;i++) {
    *file << setw(8) << setprecision(3) << showpoint << currentAtom->at(i); /* positional coordinate in Angstroem */
  }
  *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getValence(); /* occupancy */
  *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getNoValenceOrbitals(); /* temperature factor */
  *file << noshowpoint;
  *file << setw(6) << "      ";
  *file << setw(4) << "0";
  *file << setfill(' ') << setw(2) << currentAtom->getType()->getSymbol();
  *file << setw(2) << "0";

  *file << endl;
}

/**
 * Writes the neighbor information of one atom to the provided stream.
 *
 * \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;
  if (!currentAtom->ListOfBonds.empty()) {
    for(BondList::iterator currentBond = currentAtom->ListOfBonds.begin(); currentBond != currentAtom->ListOfBonds.end(); ++currentBond) {
      if (MaxNo >= MaxnumberOfNeighbors) {
        *file << "CONECT";
        *file << setw(5) << getAtomId(currentAtom->getId());
        MaxNo = 0;
      }
      *file << setw(5) << getAtomId((*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
{
  ConvertTo<size_t> toSize_t;
  ASSERT(additionalAtomData.find(atomid) != additionalAtomData.end(),
      "PdbParser::getSerial: atomid "+toString(atomid)+" not present in Map.");
  const PdbAtomInfoContainer &atomInfo = additionalAtomData.at(atomid);

  return toSize_t(atomInfo.get(PdbKey::serial));
}

/** Retrieves a value from PdbParser::atomIdMap.
 * \param atomid key
 * \return value
 */
size_t PdbParser::getAtomId(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::setAtomId(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.");
}

/** 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 line to parse as an atom
 * \param newmol molecule to add parsed atoms to
 */
void PdbParser::readAtomDataLine(std::string &line, molecule *newmol = NULL) {
  vector<string>::iterator it;
  stringstream lineStream;
  atom* newAtom = World::getInstance().createAtom();
  additionalAtomData[newAtom->getId()] = *(new PdbAtomInfoContainer);
  PdbAtomInfoContainer &atomInfo = additionalAtomData[newAtom->getId()];
  string word;
  ConvertTo<size_t> toSize_t;
  double tmp;

  lineStream << line;
  atomInfo.set(PdbKey::serial, line.substr(6,5));
  std::pair< std::set<size_t>::const_iterator, bool> Inserter =
    SerialSet.insert(toSize_t(atomInfo.get(PdbKey::serial)));
  ASSERT(Inserter.second,
      "PdbParser::readAtomDataLine() - ATOM contains entry with serial "
      +atomInfo.get(PdbKey::serial)+" already present!");
  // assign hightest+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(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);

  setAtomId(toSize_t(atomInfo.get(PdbKey::serial)), newAtom->getId());
  atomInfo.set(PdbKey::name, line.substr(12,4));
  atomInfo.set(PdbKey::altloc, line.substr(16,1));
  atomInfo.set(PdbKey::resName, line.substr(17,3));
  atomInfo.set(PdbKey::chainID, line.substr(21,1));
  atomInfo.set(PdbKey::resSeq, line.substr(22,4));
  atomInfo.set(PdbKey::iCode, line.substr(26,1));
  PdbAtomInfoContainer::ScanKey(tmp, line.substr(30,8));
  newAtom->set(0, tmp);
  PdbAtomInfoContainer::ScanKey(tmp, line.substr(38,8));
  newAtom->set(1, tmp);
  PdbAtomInfoContainer::ScanKey(tmp, line.substr(46,8));
  newAtom->set(2, tmp);
  atomInfo.set(PdbKey::occupancy, line.substr(54,6));
  atomInfo.set(PdbKey::tempFactor, line.substr(60,6));
  atomInfo.set(PdbKey::charge, line.substr(78,2));
  PdbAtomInfoContainer::ScanKey(word, line.substr(76,2));
  newAtom->setType(World::getInstance().getPeriode()->FindElement(word));

  if (newmol != NULL)
    newmol->AddAtom(newAtom);

//  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) << "\tserial is " << atomInfo.get(PdbKey::serial) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tname is " << atomInfo.get(PdbKey::name) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\taltloc is " << atomInfo.get(PdbKey::altloc) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tresName is " << atomInfo.get(PdbKey::resName) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tchainID is " << atomInfo.get(PdbKey::chainID) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tresSeq is " << atomInfo.get(PdbKey::resSeq) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tiCode is " << atomInfo.get(PdbKey::iCode) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tx is " << newAtom->getPosition() << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\toccupancy is " << atomInfo.get(PdbKey::occupancy) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\ttempFactor is " << atomInfo.get(PdbKey::tempFactor) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\telement is '" << *(newAtom->getType()) << "'" << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tcharge is " << atomInfo.get(PdbKey::charge) << std::endl);
}

/**
 * Reads neighbor information for one atom from the input.
 *
 * \param line to parse as an atom
 */
void PdbParser::readNeighbors(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(getAtomId(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(getAtomId(*iter)));
    _atom->addBond(_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(PdbKey::serial) == OtheratomInfo.get(PdbKey::serial));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in serials!" << std::endl);
      status = status && (atomInfo.get(PdbKey::name) == OtheratomInfo.get(PdbKey::name));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in names!" << std::endl);
      status = status && (atomInfo.get(PdbKey::altloc) == OtheratomInfo.get(PdbKey::altloc));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in altlocs!" << std::endl);
      status = status && (atomInfo.get(PdbKey::resName) == OtheratomInfo.get(PdbKey::resName));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resNames!" << std::endl);
      status = status && (atomInfo.get(PdbKey::chainID) == OtheratomInfo.get(PdbKey::chainID));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in chainIDs!" << std::endl);
      status = status && (atomInfo.get(PdbKey::resSeq) == OtheratomInfo.get(PdbKey::resSeq));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resSeqs!" << std::endl);
      status = status && (atomInfo.get(PdbKey::iCode) == OtheratomInfo.get(PdbKey::iCode));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in iCodes!" << std::endl);
      status = status && (atomInfo.get(PdbKey::occupancy) == OtheratomInfo.get(PdbKey::occupancy));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in occupancies!" << std::endl);
      status = status && (atomInfo.get(PdbKey::tempFactor) == OtheratomInfo.get(PdbKey::tempFactor));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in tempFactors!" << std::endl);
      status = status && (atomInfo.get(PdbKey::charge) == OtheratomInfo.get(PdbKey::charge));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in charges!" << std::endl);
    }
  }

  return status;
}