source: src/Parser/XyzParser.cpp@ bbff92

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

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

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

#include "CodePatterns/MemDebug.hpp"

#include <vector>

#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"

#include "XyzParser.hpp"

#include "atom.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "World.hpp"

using namespace std;

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

/**
 * Constructor.
 */
FormatParser< xyz >::FormatParser() :
  FormatParser_common(NULL),
  comment("")
{}

/**
 * Destructor.
 */
FormatParser< xyz >::~FormatParser() 
{}

/**
 * Loads an XYZ file into the World.
 *
 * \param XYZ file
 */
void FormatParser< xyz >::load(istream* file)
{
  atom* newAtom = NULL;
  molecule* newmol = NULL;
  int numberOfAtoms;
  char commentBuffer[512], type[3];
  string number;
  std::vector<atom *> AddedAtoms;

  // create the 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);

  // the first line tells number of atoms,
  // where we need this construct due to skipping of empty lines below
  *file >> number >> ws;
  unsigned int step = 0;
  while (file->good()) {
    std::stringstream numberstream(number);
    numberstream >> numberOfAtoms;
    if (step == 0)
      AddedAtoms.resize(numberOfAtoms);
    // the second line is always a comment
    file->getline(commentBuffer, 512);
    comment = commentBuffer;

    for (int i = 0; i < numberOfAtoms; i++) {
      // parse type and position
      *file >> type;
      Vector tempVector;
      for (int j=0;j<NDIM;j++) {
        *file >> tempVector[j];
      }
      LOG(4, "INFO: Parsed type as " << type << " and position at "
          << tempVector << " for time step " << step);

      if (step == 0) {
        newAtom = World::getInstance().createAtom();
        newAtom->setType(World::getInstance().getPeriode()->FindElement(type));
        newmol->AddAtom(newAtom);
        AddedAtoms[i] = newAtom;
        LOG(5, "INFO: Created new atom, setting " << i << " th component of AddedAtoms.");
      } else {
        newAtom = AddedAtoms[i];
        LOG(5, "INFO: Using present atom " << *newAtom << " from " << i << " th component of AddedAtoms.");
        ASSERT(newAtom->getType() == World::getInstance().getPeriode()->FindElement(type),
            "FormatParser< xyz >::load() - atom has different type "+newAtom->getType()->getSymbol()+" than parsed now "+type+", mixed up order?");
      }
      newAtom->setPositionAtStep(step, tempVector);
    }
    getline (*file, number);  // eat away rest of last line
    // skip empty lines
    unsigned int counter = 0;
    while (file->good()) {
      getline (*file, number);
      LOG(4, "INFO: Skipped line is '" << number << "'");
      counter++;
      if (!number.empty())
        break;
    }
    LOG(3, "INFO: I skipped "+toString(counter)+" empty lines.");
    ++step;
  }
  BOOST_FOREACH(atom *_atom, World::getInstance().getAllAtoms())
    LOG(3, "INFO: Atom " << _atom->getName() << " " << *dynamic_cast<AtomInfo *>(_atom) << ".");

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

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

  // get max and min trajectories
  unsigned int min_trajectories = 1;
  unsigned int max_trajectories = 1;
  for (std::vector<atom *>::const_iterator iter = atoms.begin();
      iter != atoms.end();
      ++iter) {
    if (max_trajectories < (*iter)->getTrajectorySize())
      max_trajectories = (*iter)->getTrajectorySize();
    if ((min_trajectories > (*iter)->getTrajectorySize())
        && (max_trajectories > (*iter)->getTrajectorySize()))
      min_trajectories = (*iter)->getTrajectorySize();
  }
  ASSERT((min_trajectories == max_trajectories) || (min_trajectories == 1),
      "FormatParser< xyz >::save() - not all atoms have same number of trajectories.");
  LOG(2, "INFO: There are " << max_trajectories << " to save.");

  for (unsigned int step = 0; step < max_trajectories; ++step) {
    if (step != 0)
      *file << "\n";
    //if (comment == "") {
      time_t now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
      comment = "Created by molecuilder on ";
      // 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)
        comment += time.substr(0,pos);
      else
        comment += time;
      comment += ", time step "+toString(step);
    //}
    *file << atoms.size() << endl << "\t" << comment << endl;

    for(vector<atom*>::const_iterator it = atoms.begin(); it != atoms.end(); it++) {
      *file << noshowpoint << (*it)->getType()->getSymbol();
      *file << "\t" << (*it)->atStep(0, step);
      *file << "\t" << (*it)->atStep(1, step);
      *file << "\t" << (*it)->atStep(2, step);
      *file << endl;
    }
  }
}