source: src/Parser/PdbParser.cpp@ a42054

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

#include <map>
#include <vector>

#include <iostream>
#include <iomanip>

using namespace std;

/**
 * Constructor.
 */
PdbParser::PdbParser() {
  knownKeys[" "] = PdbKey::noKey; // with this we can detect invalid keys
  knownKeys["x"] = PdbKey::x;
  knownKeys["Id"] = PdbKey::Id;
  knownKeys["Type"] = PdbKey::Type;
  knownKeys["extType"] = PdbKey::extType;
  knownKeys["name"] = PdbKey::name;
  knownKeys["resName"] = PdbKey::resName;
  knownKeys["chainID"] = PdbKey::chainID;
  knownKeys["resSeq"] = PdbKey::resSeq;
  knownKeys["occupancy"] = PdbKey::occupancy;
  knownKeys["tempFactor"] = PdbKey::tempFactor;
  knownKeys["segID"] = PdbKey::segID;
  knownKeys["charge"] = PdbKey::charge;
}

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

/**
 * Loads atoms from a tremolo-formatted file.
 *
 * \param tremolo file
 */
void PdbParser::load(istream* file) {
//  string line;
//  string::size_type location;
//
//  usedFields.clear();
//  while (file->good()) {
//    std::getline(*file, line, '\n');
//    if (usedFields.empty()) {
//      location = line.find("ATOMDATA", 0);
//      if (location != string::npos) {
//       parseAtomDataKeysLine(line, location + 8);
//      }
//    }
//    if (line.length() > 0 && line.at(0) != '#') {
//      readAtomDataLine(line);
//    }
//  }
//
//  processNeighborInformation();
//  adaptImprData();
//  adaptTorsion();
}

/**
 * Saves the World's current state into as a tremolo 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;
  }

  {
    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 = 0 is reserved
        MolNo = 0;
      } 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 tremolo-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) {
  if (!currentAtom->ListOfBonds.empty()) {
    *file << "CONECT";
    *file << setw(5) << getAtomId(currentAtom->getId());
    int MaxNo = 0;
    for(BondList::iterator currentBond = currentAtom->ListOfBonds.begin(); currentBond != currentAtom->ListOfBonds.end(); ++currentBond) {
      if (MaxNo < MaxnumberOfNeighbors) {
        *file << setw(5) << getAtomId((*currentBond)->GetOtherAtom(currentAtom)->getId());
      }
      MaxNo++;
    }
    *file << endl;
  }
}


/** Retrieves a value from PdbParser::atomIdMap.
 * \param atomid key
 * \return value
 */
int PdbParser::getAtomId(int 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(int localatomid, int atomid)
{
  pair<std::map<int,int>::iterator, bool > inserter;
  inserter = atomIdMap.insert( pair<int, int>(localatomid, atomid) );
  ASSERT(inserter.second, "PdbParser::setAtomId: atomId already present in Map.");
}

/**
 * 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
 */
void PdbParser::readAtomDataLine(string line) {
//  vector<string>::iterator it;
//  stringstream lineStream;
//  atom* newAtom = World::getInstance().createAtom();
//  PdbAtomInfoContainer *atomInfo = NULL;
//  additionalAtomData[newAtom->getId()] = *(new PdbAtomInfoContainer);
//  atomInfo = &additionalAtomData[newAtom->getId()];
//  PdbKey::atomDataKey currentField;
//  string word;
//  int oldId;
//  double tmp;
//
//  lineStream << line;
//  for (it = usedFields.begin(); it < usedFields.end(); it++) {
//    currentField = knownKeys[it->substr(0, it->find("="))];
//    switch (currentField) {
//      case PdbKey::x :
//        // for the moment, assume there are always three dimensions
//        for (int i=0;i<NDIM;i++) {
//          lineStream >> tmp;
//          newAtom->set(i, tmp);
//        }
//        break;
//      case PdbKey::u :
//        // for the moment, assume there are always three dimensions
//        lineStream >> newAtom->AtomicVelocity[0];
//        lineStream >> newAtom->AtomicVelocity[1];
//        lineStream >> newAtom->AtomicVelocity[2];
//        break;
//      case PdbKey::Type :
//        char type[3];
//        lineStream >> type;
//        newAtom->setType(World::getInstance().getPeriode()->FindElement(type));
//        ASSERT(newAtom->getType(), "Type was not set for this atom");
//        break;
//      case PdbKey::Id :
//        lineStream >> oldId;
//        atomIdMap[oldId] = newAtom->getId();
//        break;
//      case PdbKey::neighbors :
//        readNeighbors(&lineStream,
//            atoi(it->substr(it->find("=") + 1, 1).c_str()), newAtom->getId());
//        break;
//      default :
//        lineStream >> word;
//        atomInfo->set(currentField, word);
//        break;
//    }
//  }
}

/**
 * Reads neighbor information for one atom from the input.
 *
 * \param stream where to read the information from
 * \param number of neighbors to read
 * \param world id of the atom the information belongs to
 */
void PdbParser::readNeighbors(stringstream* line, int numberOfNeighbors, 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) {
//      additionalAtomData[atomId].neighbors.push_back(neighborId);
//    }
//  }
}

/**
 * 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.
 */
void PdbParser::processNeighborInformation() {
//  if (!isUsedField("neighbors")) {
//    return;
//  }
//
//  for(map<int, PdbAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
//    currentInfo != additionalAtomData.end(); currentInfo++
//  ) {
//    for(vector<int>::iterator neighbor = currentInfo->second.neighbors.begin();
//      neighbor != currentInfo->second.neighbors.end(); neighbor++
//    ) {
//      World::getInstance().getAtom(AtomById(currentInfo->first))
//          ->addBond(World::getInstance().getAtom(AtomById(atomIdMap[*neighbor])));
//    }
//  }
}

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


PdbAtomInfoContainer::PdbAtomInfoContainer() :
  name("-"),
  resName("-"),
  chainID("0"),
  resSeq("0"),
  occupancy("0"),
  tempFactor("0"),
  segID("0"),
  charge("0")
{}

void PdbAtomInfoContainer::set(PdbKey::PdbDataKey key, string value) {
  switch (key) {
    case PdbKey::extType :
      extType = value;
      break;
    case PdbKey::name :
      name = value;
      break;
    case PdbKey::resName :
      resName = value;
      break;
    case PdbKey::chainID :
      chainID = value;
      break;
    case PdbKey::resSeq :
      resSeq = value;
      break;
    case PdbKey::occupancy :
      occupancy = value;
      break;
    case PdbKey::tempFactor :
      tempFactor = value;
      break;
    case PdbKey::segID :
      segID = value;
      break;
    case PdbKey::charge :
      charge = value;
      break;
    default :
      cout << "Unknown key: " << key << ", value: " << value << endl;
      break;
  }
}

string PdbAtomInfoContainer::get(PdbKey::PdbDataKey key) {
  switch (key) {
    case PdbKey::extType :
      return extType;
    case PdbKey::name :
      return name;
    case PdbKey::resName :
      return resName;
    case PdbKey::chainID :
      return chainID;
    case PdbKey::resSeq :
      return resSeq;
    case PdbKey::occupancy :
      return occupancy;
    case PdbKey::tempFactor :
      return tempFactor;
    case PdbKey::segID :
      return segID;
    case PdbKey::charge :
      return charge;
    default :
      cout << "Unknown key: " << key << endl;
      return "";
  }
}