source: src/Parser/XyzParser.cpp@ 1e6591

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

/*
 * 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 <limits>
#include <vector>

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

#include "XyzParser.hpp"

#include "Atom/atom.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "molecule.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;
  string elementtype;
  string streambuffer;
  std::vector<atom *> AddedAtoms;

  // create the molecule
  newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;

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

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

      if (step == 0) {
        newAtom = World::getInstance().createAtom();
        newAtom->setType(World::getInstance().getPeriode()->FindElement(elementtype.c_str()));
        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(elementtype.c_str()),
            "FormatParser< xyz >::load() - atom has different type "+newAtom->getType()->getSymbol()+" than parsed now "+elementtype+", mixed up order?");
      }
      newAtom->setPositionAtStep(step, tempVector);
    }
    getline (*file, streambuffer);  // eat away rest of last line
    // skip empty lines
    unsigned int counter = 0;
    while (file->good()) {
      getline (*file, streambuffer);
      LOG(4, "INFO: Skipped line is '" << streambuffer << "'");
      counter++;
      if (!streambuffer.empty())
        break;
    }
    LOG(3, "INFO: I skipped "+toString(counter)+" empty lines.");
    ++step;
  }
  BOOST_FOREACH(const atom *_atom, const_cast<const World &>(World::getInstance()).getAllAtoms())
    LOG(3, "INFO: Atom " << _atom->getName() << " " << *dynamic_cast<const AtomInfo *>(_atom) << ".");

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

const std::string FormatParser< xyz >::printCoordinate(const double value)
{
  std::stringstream position;
  if (fabs(value) > 1000) // enforce precision for large components
    position << std::fixed << std::setprecision(3) << value; 
  else
    position << value; 
  return position.str();
}

/**
 * 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<const atom *> &atoms) {
  LOG(2, "DEBUG: Saving changes to xyz.");

  // get max and min trajectories
  std::pair<size_t, size_t> minmax_trajectories =
      getMinMaxTrajectories(atoms);
  LOG(2, "INFO: There are " << minmax_trajectories.second << " steps to save.");

  // always store at least one step
  for (size_t step = 0; (step < minmax_trajectories.second) || (step == 0); ++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<const atom*>::const_iterator it = atoms.begin(); it != atoms.end(); it++) {
      *file << (*it)->getType()->getSymbol();
      *file << "\t" << printCoordinate((*it)->atStep(0, step));
      *file << "\t" << printCoordinate((*it)->atStep(1, step));
      *file << "\t" << printCoordinate((*it)->atStep(2, step));
      *file << endl;
    }
  }
}