source: src/Parser/TremoloParser.cpp@ 56f73b

/*
 * 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.
 */

/*
 * TremoloParser.cpp
 *
 *  Created on: Mar 2, 2010
 *      Author: metzler
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "TremoloParser.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.
 */
TremoloParser::TremoloParser() {
  knownKeys["x"] = TremoloKey::x;
  knownKeys["u"] = TremoloKey::u;
  knownKeys["F"] = TremoloKey::F;
  knownKeys["stress"] = TremoloKey::stress;
  knownKeys["Id"] = TremoloKey::Id;
  knownKeys["neighbors"] = TremoloKey::neighbors;
  knownKeys["imprData"] = TremoloKey::imprData;
  knownKeys["GroupMeasureTypeNo"] = TremoloKey::GroupMeasureTypeNo;
  knownKeys["Type"] = TremoloKey::Type;
  knownKeys["extType"] = TremoloKey::extType;
  knownKeys["name"] = TremoloKey::name;
  knownKeys["resName"] = TremoloKey::resName;
  knownKeys["chainID"] = TremoloKey::chainID;
  knownKeys["resSeq"] = TremoloKey::resSeq;
  knownKeys["occupancy"] = TremoloKey::occupancy;
  knownKeys["tempFactor"] = TremoloKey::tempFactor;
  knownKeys["segID"] = TremoloKey::segID;
  knownKeys["Charge"] = TremoloKey::Charge;
  knownKeys["charge"] = TremoloKey::charge;
  knownKeys["GrpTypeNo"] = TremoloKey::GrpTypeNo;
  knownKeys["torsion"] = TremoloKey::torsion;

  // default behavior: use all possible keys on output
  for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin(); iter != knownKeys.end(); ++iter)
    usedFields.push_back(iter->first);

  // and noKey afterwards(!) such that it is not used in usedFields
  knownKeys[" "] = TremoloKey::noKey; // with this we can detect invalid keys
}

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

/**
 * Loads atoms from a tremolo-formatted file.
 *
 * \param tremolo file
 */
void TremoloParser::load(istream* file) {
  string line;
  string::size_type location;

  // reset atomIdMap, for we now get new serials
  atomIdMap.clear();
  usedFields.clear();

  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
  World::getInstance().getMolecules()->insert(newmol);
  while (file->good()) {
    std::getline(*file, line, '\n');
    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, newmol);
    }
  }

  processNeighborInformation();
  adaptImprData();
  adaptTorsion();
}

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

  vector<atom*>::const_iterator atomIt;
  vector<string>::iterator it;

  *file << "# ATOMDATA";
  for (it=usedFields.begin(); it < usedFields.end(); it++) {
    *file << "\t" << *it;
  }
  *file << endl;
  for (atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
    saveLine(file, *atomIt);
  }
}

/**
 * Sets the keys for which data should be written to the stream when save is
 * called.
 *
 * \param string of field names with the same syntax as for an ATOMDATA line
 *        but without the prexix "ATOMDATA"
 */
void TremoloParser::setFieldsForSave(std::string atomDataLine) {
  parseAtomDataKeysLine(atomDataLine, 0);
}


/**
 * Writes one line of tremolo-formatted data to the provided stream.
 *
 * \param stream where to write the line to
 * \param reference to the atom of which information should be written
 */
void TremoloParser::saveLine(ostream* file, atom* currentAtom) {
  vector<string>::iterator it;
  TremoloKey::atomDataKey currentField;

  for (it = usedFields.begin(); it != usedFields.end(); it++) {
    currentField = knownKeys[it->substr(0, it->find("="))];
    switch (currentField) {
      case TremoloKey::x :
        // for the moment, assume there are always three dimensions
        *file << currentAtom->at(0) << "\t";
        *file << currentAtom->at(1) << "\t";
        *file << currentAtom->at(2) << "\t";
        break;
      case TremoloKey::u :
        // for the moment, assume there are always three dimensions
        *file << currentAtom->AtomicVelocity[0] << "\t";
        *file << currentAtom->AtomicVelocity[1] << "\t";
        *file << currentAtom->AtomicVelocity[2] << "\t";
        break;
      case TremoloKey::Type :
        *file << currentAtom->getType()->getSymbol() << "\t";
        break;
      case TremoloKey::Id :
        *file << currentAtom->getId()+1 << "\t";
        break;
      case TremoloKey::neighbors :
        writeNeighbors(file, atoi(it->substr(it->find("=") + 1, 1).c_str()), currentAtom);
        break;
      case TremoloKey::resSeq :
        if (additionalAtomData.find(currentAtom->getId()) != additionalAtomData.end()) {
          *file << additionalAtomData[currentAtom->getId()].get(currentField);
        } else if (currentAtom->getMolecule() != NULL) {
          *file << setw(4) << currentAtom->getMolecule()->getId()+1;
        } else {
          *file << defaultAdditionalData.get(currentField);
        }
        *file << "\t";
      break;
      default :
        if (additionalAtomData.find(currentAtom->getId()) != additionalAtomData.end()) {
          *file << additionalAtomData[currentAtom->getId()].get(currentField);
        } else if (additionalAtomData.find(currentAtom->GetTrueFather()->getId()) != additionalAtomData.end()) {
          *file << additionalAtomData[currentAtom->GetTrueFather()->getId()].get(currentField);
        } else {
          *file << defaultAdditionalData.get(currentField);
        }
        *file << "\t";
        break;
    }
  }

  *file << endl;
}

/**
 * Writes the neighbor information of one atom to the provided stream.
 *
 * \param stream where to write neighbor information to
 * \param number of neighbors
 * \param reference to the atom of which to take the neighbor information
 */
void TremoloParser::writeNeighbors(ostream* file, int numberOfNeighbors, atom* currentAtom) {
  BondList::iterator currentBond = currentAtom->ListOfBonds.begin();
  for (int i = 0; i < numberOfNeighbors; i++) {
    *file << (currentBond != currentAtom->ListOfBonds.end()
        ? (*currentBond)->GetOtherAtom(currentAtom)->getId()+1 : 0) << "\t";
    if (currentBond != currentAtom->ListOfBonds.end())
      ++currentBond;
  }
}

/**
 * Stores keys from the ATOMDATA line.
 *
 * \param line to parse the keys from
 * \param with which offset the keys begin within the line
 */
void TremoloParser::parseAtomDataKeysLine(string line, int offset) {
  string keyword;
  stringstream lineStream;

  lineStream << line.substr(offset);
  usedFields.clear();
  while (lineStream.good()) {
    lineStream >> keyword;
    if (knownKeys[keyword.substr(0, keyword.find("="))] == TremoloKey::noKey) {
      // TODO: throw exception about unknown key
      cout << "Unknown key: " << keyword << " is not part of the tremolo format specification." << endl;
      break;
    }
    usedFields.push_back(keyword);
  }
}

/**
 * Reads one data line of a tremolo file and interprets it according to the keys
 * obtained from the ATOMDATA line.
 *
 * \param line to parse as an atom
 * \param *newmol molecule to add atom to
 */
void TremoloParser::readAtomDataLine(string line, molecule *newmol = NULL) {
  vector<string>::iterator it;
  stringstream lineStream;
  atom* newAtom = World::getInstance().createAtom();
  TremoloAtomInfoContainer *atomInfo = NULL;
  additionalAtomData[newAtom->getId()] = TremoloAtomInfoContainer(); // fill with default values
  atomInfo = &additionalAtomData[newAtom->getId()];
  TremoloKey::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 TremoloKey::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 TremoloKey::u :
        // for the moment, assume there are always three dimensions
        lineStream >> newAtom->AtomicVelocity[0];
        lineStream >> newAtom->AtomicVelocity[1];
        lineStream >> newAtom->AtomicVelocity[2];
        break;
      case TremoloKey::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 TremoloKey::Id :
        lineStream >> oldId;
        atomIdMap[oldId] = newAtom->getId();
        break;
      case TremoloKey::neighbors :
        readNeighbors(&lineStream,
            atoi(it->substr(it->find("=") + 1, 1).c_str()), newAtom->getId());
        break;
      default :
        lineStream >> word;
        atomInfo->set(currentField, word);
        break;
    }
  }
  if (newmol != NULL)
    newmol->AddAtom(newAtom);
}

/**
 * Reads neighbor information for one atom from the input.
 *
 * \param line stream where to read the information from
 * \param numberOfNeighbors number of neighbors to read
 * \param atomid world id of the atom the information belongs to
 */
void TremoloParser::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) {
//      std::cout << "Atom with global id " << atomId << " has neighbour with serial " << neighborId << std::endl;
      additionalAtomData[atomId].neighbors.push_back(neighborId);
    }
  }
}

/**
 * Checks whether the provided name is within the list of used fields.
 *
 * \param field name to check
 *
 * \return true if the field name is used
 */
bool TremoloParser::isUsedField(string fieldName) {
  bool fieldNameExists = false;
  for (vector<string>::iterator usedField = usedFields.begin(); usedField != usedFields.end(); usedField++) {
    if (usedField->substr(0, usedField->find("=")) == fieldName)
      fieldNameExists = true;
  }

  return fieldNameExists;
}


/**
 * Adds the collected neighbor information to the atoms in the world. The atoms
 * are found by their current ID and mapped to the corresponding atoms with the
 * Id found in the parsed file.
 */
void TremoloParser::processNeighborInformation() {
  if (!isUsedField("neighbors")) {
    return;
  }

  for(map<int, TremoloAtomInfoContainer>::iterator currentInfo = additionalAtomData.begin();
    currentInfo != additionalAtomData.end(); currentInfo++
  ) {
    if (!currentInfo->second.neighbors_processed) {
      for(vector<int>::iterator neighbor = currentInfo->second.neighbors.begin();
        neighbor != currentInfo->second.neighbors.end(); neighbor++
      ) {
//        std::cout << "Creating bond between ("
//            << currentInfo->first
//            << ") and ("
//            << atomIdMap[*neighbor] << "|" << *neighbor << ")" << std::endl;
        World::getInstance().getAtom(AtomById(currentInfo->first))
            ->addBond(World::getInstance().getAtom(AtomById(atomIdMap[*neighbor])));
      }
      currentInfo->second.neighbors_processed = true;
    }
  }
}

/**
 * Replaces atom IDs read from the file by the corresponding world IDs. All IDs
 * IDs of the input string will be replaced; expected separating characters are
 * "-" and ",".
 *
 * \param string in which atom IDs should be adapted
 *
 * \return input string with modified atom IDs
 */
string TremoloParser::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();
}

/**
 * Corrects the atom IDs in each imprData entry to the corresponding world IDs
 * as they might differ from the originally read IDs.
 */
void TremoloParser::adaptImprData() {
  if (!isUsedField("imprData")) {
    return;
  }

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

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

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