/*
 * OutputEnergies.hpp
 *
 *  Created on: Feb 23, 2011
 *      Author: heber
 */

#ifndef OUTPUTTEMPERATURE_HPP_
#define OUTPUTTEMPERATURE_HPP_

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

#include "Helpers/defs.hpp"

template <class T>
class OutputEnergies
{
public:
  OutputEnergies(AtomSetMixin<T> &_atoms) :
    atoms(_atoms)
  {}
  ~OutputEnergies()
  {}

  /** Stores the temperature evaluated from velocities in molecule::Trajectories.
   * We simply use the formula equivaleting temperature and kinetic energy:
   * \f$k_B T = \sum_i m_i v_i^2\f$
   * \param *output output stream of temperature file
   * \param startstep first MD step in molecule::Trajectories
   * \param endstep last plus one MD step in molecule::Trajectories
   * \return file written (true), failure on writing file (false)
   */
  bool operator()(ofstream * const output, int startstep, int endstep)
  {
    double temperature;
    Vector force, abs_force;
    Vector momentum, abs_momentum;
    // test stream
    if (output == NULL)
      return false;
    else
      *output << "# Step"
        << "\tTemperature [K]"
        << "\tTemperature [a.u.]"
        << "\tMomentum"
        << "\tAbolute Momentum"
        << "\tForce"
        << "\tAbsolute Force"
        << endl;
    for (int step=startstep;step < endstep; step++) { // loop over all time steps
      temperature = atoms.totalTemperatureAtStep(step);
      momentum = atoms.totalMomentumAtStep(step);
      force = atoms.totalForceAtStep(step);
      abs_momentum = atoms.totalAbsoluteMomentumAtStep(step);
      abs_force = atoms.totalAbsoluteForceAtStep(step);
      *output << step
          << "\t" << temperature*AtomicMassUnitsAngstroemOverAtomictimeSquaredToHt*AtomicEnergyToKelvin
          << "\t" << temperature*AtomicMassUnitsAngstroemOverAtomictimeSquaredToHt
          << "\t" << momentum.Norm()*AtomicMassUnitsAngstroemOverAtomictimeToAtomicMomentum
          << "\t" << abs_momentum.Norm()*AtomicMassUnitsAngstroemOverAtomictimeToAtomicMomentum
          << "\t" << force.Norm()
          << "\t" << abs_force.Norm()
          << endl;
    }
    return true;
  };

private:
  AtomSetMixin<T> atoms;
};

#endif /* OUTPUTTEMPERATURE_HPP_ */