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

#ifndef ATOMICFORCEMANIPULATOR_HPP_
#define ATOMICFORCEMANIPULATOR_HPP_

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

#include "Atom/atom.hpp"
#include "Atom/AtomSet.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Helpers/helpers.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "Thermostats/ThermoStatContainer.hpp"
#include "Thermostats/Thermostat.hpp"
#include "World.hpp"

template <class T>
class AtomicForceManipulator
{
public:
  /** Constructor of class AtomicForceManipulator.
   *
   * \param _atoms set of atoms to integrate
   * \param _Deltat time step width in atomic units
   * \param _IsAngstroem whether length units are in angstroem or bohr radii
   */
  AtomicForceManipulator(AtomSetMixin<T> &_atoms, double _Deltat, bool _IsAngstroem) :
    Deltat(_Deltat),
    IsAngstroem(_IsAngstroem),
    atoms(_atoms)
  {}
  /** Destructor of class AtomicForceManipulator.
   *
   */
  ~AtomicForceManipulator()
  {}

  /** Parses nuclear forces from file.
   * Forces are stored in the time step \a TimeStep in the atomicForces in \a atoms.
   * \param *file filename
   * \param TimeStep time step to parse forces file into
   * \return true - file parsed, false - file not found or imparsable
   */
  bool parseForcesFile(const char *file, const int TimeStep)
  {
    Info FunctionInfo(__func__);
    ForceMatrix Force;

    // parse file into ForceMatrix
    std::ifstream input(file);
    if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
      ELOG(0, "Could not parse Force Matrix file " << file << ".");
      return false;
    }
    input.close();
    if (Force.RowCounter[0] != (int)atoms.size()) {
      ELOG(0, "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << atoms.size() << ".");
      return false;
    }

    addForceMatrixToAtomicForce(Force, TimeStep, 1);
    return true;
  }

protected:
  void addForceMatrixToAtomicForce(const ForceMatrix &Force, const int &TimeStep, const int offset)
  {
    // place forces from matrix into atoms
    Vector tempVector;
    size_t i=0;
    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter,++i) {
      for(size_t d=0;d<NDIM;d++) {
        tempVector[d] = Force.Matrix[0][i][d+offset]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
      }
      tempVector += (*iter)->getAtomicForceAtStep(TimeStep);
      (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
    }
  }

  void correctForceMatrixForFixedCenterOfMass(const size_t offset, const int &TimeStep) {
    Vector ForceVector;
    // correct Forces
    //std::cout << "Force before correction, " << Force << std::endl;
    ForceVector.Zero();
    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
      ForceVector += (*iter)->getAtomicForceAtStep(TimeStep);
    }
    ForceVector.Scale(1./(double)atoms.size());
    //std::cout << "Force before second correction, " << Force << std::endl;
    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
      const Vector tempVector = (*iter)->getAtomicForceAtStep(TimeStep) - ForceVector;
      (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
    }
    LOG(3, "INFO: forces corrected by " << ForceVector << " each.");
  }

  void correctVelocitiesForFixedCenterOfMass(const int &TimeStep) {
    Vector Velocity;
    double IonMass;
    // correct velocities (rather momenta) so that center of mass remains motionless
    Velocity = atoms.totalMomentumAtStep(TimeStep);
    IonMass = atoms.totalMass();

    // correct velocities (rather momenta) so that center of mass remains motionless
    Velocity *= 1./IonMass;
    atoms.addVelocityAtStep(-1.*Velocity,TimeStep);

    LOG(3, "INFO: Velocities corrected by " << Velocity << " each.");
  }

  void performThermostatControl(const int &TimeStep) {
    double ActualTemp;

    // calculate current temperature
    ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
    LOG(3, "INFO: Current temperature is " << ActualTemp);

    // rescale to desired value
    double ekin = World::getInstance().getThermostats()->getActive()->scaleAtoms(TimeStep,ActualTemp,atoms);
    ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
    LOG(3, "INFO: New temperature after thermostat is " << ActualTemp);
    LOG(1, "Kinetic energy is " << ekin << ".");
  }

protected:
  double Deltat;
  bool IsAngstroem;
  AtomSetMixin<T> atoms;
};

#endif /* ATOMICFORCEMANIPULATOR_HPP_ */