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

/*
 * molecule_dynamics.cpp
 *
 *  Created on: Oct 5, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "World.hpp"
#include "atom.hpp"
#include "config.hpp"
#include "Dynamics/MinimiseConstrainedPotential.hpp"
//#include "element.hpp"
#include "CodePatterns/Info.hpp"
//#include "CodePatterns/Verbose.hpp"
//#include "CodePatterns/Log.hpp"
#include "molecule.hpp"
#include "parser.hpp"
//#include "LinearAlgebra/Plane.hpp"
#include "ThermoStatContainer.hpp"
#include "Thermostats/Berendsen.hpp"

#include "CodePatterns/enumeration.hpp"

/************************************* Functions for class molecule *********************************/


/** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
 * Note, step number is config::MaxOuterStep
 * \param *out output stream for debugging
 * \param startstep stating initial configuration in molecule::Trajectories
 * \param endstep stating final configuration in molecule::Trajectories
 * \param &prefix path and prefix
 * \param &config configuration structure
 * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential()
 * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
 */
bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, std::string prefix, config &configuration, bool MapByIdentity)
{
  // TODO: rewrite permutationMaps using enumeration objects
  molecule *mol = NULL;
  bool status = true;
  int MaxSteps = configuration.MaxOuterStep;
  MoleculeListClass *MoleculePerStep = new MoleculeListClass(World::getPointer());
  // Get the Permutation Map by MinimiseConstrainedPotential
  atom **PermutationMap = NULL;
  atom *Sprinter = NULL;
  if (!MapByIdentity) {
    MinimiseConstrainedPotential Minimiser(atoms);
    Minimiser(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
  } else {
    // TODO: construction of enumeration goes here
    PermutationMap = new atom *[getAtomCount()];
    for(internal_iterator iter = atoms.begin(); iter != atoms.end();++iter){
      PermutationMap[(*iter)->getNr()] = (*iter);
    }
  }

  // check whether we have sufficient space in Trajectories for each atom
  LOG(1, "STATUS: Extending each trajectory size to " << MaxSteps+1 << ".");
  for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::ResizeTrajectory),MaxSteps+1));
  // push endstep to last one
  for_each(atoms.begin(),atoms.end(),boost::bind(&atom::CopyStepOnStep,_1,MaxSteps,endstep));
  endstep = MaxSteps;

  // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
  DoLog(1) && (Log() << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl);
  for (int step = 0; step <= MaxSteps; step++) {
    mol = World::getInstance().createMolecule();
    MoleculePerStep->insert(mol);
    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
      // add to molecule list
      Sprinter = mol->AddCopyAtom((*iter));
      // add to Trajectories
      Vector temp = (*iter)->getPositionAtStep(startstep) + (PermutationMap[(*iter)->getNr()]->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep))*((double)step/(double)MaxSteps);
      Sprinter->setPosition(temp);
      (*iter)->setAtomicVelocityAtStep(step, zeroVec);
      (*iter)->setAtomicForceAtStep(step, zeroVec);
      //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    }
  }
  MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit

  // store the list to single step files
  int *SortIndex = new int[getAtomCount()];
  for (int i=getAtomCount(); i--; )
    SortIndex[i] = i;

  status = MoleculePerStep->OutputConfigForListOfFragments(prefix, SortIndex);
  delete[](SortIndex);

  // free and return
  delete[](PermutationMap);
  delete(MoleculePerStep);
  return status;
};

/** Parses nuclear forces from file and performs Verlet integration.
 * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
 * have to transform them).
 * This adds a new MD step to the config file.
 * \param *file filename
 * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
 * \param offset offset in matrix file to the first force component
 * \return true - file found and parsed, false - file not found or imparsable
 * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
 */
bool molecule::VerletForceIntegration(char *file, config &configuration, const size_t offset)
{
  Info FunctionInfo(__func__);
  string token;
  stringstream item;
  double IonMass, ConstrainedPotentialEnergy, ActualTemp;
  Vector Velocity;
  ForceMatrix Force;

  const int AtomCount = getAtomCount();
  // parse file into ForceMatrix
  std::ifstream input(file);
  if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
    DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
    performCriticalExit();
    return false;
  }
  input.close();
  if (Force.RowCounter[0] != AtomCount) {
    DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
    performCriticalExit();
    return false;
  }
  // correct Forces
  Velocity.Zero();
  for(int i=0;i<AtomCount;i++)
    for(int d=0;d<NDIM;d++) {
      Velocity[d] += Force.Matrix[0][i][d+offset];
    }
  for(int i=0;i<AtomCount;i++)
    for(int d=0;d<NDIM;d++) {
      Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount);
    }
  // solve a constrained potential if we are meant to
  if (configuration.DoConstrainedMD) {
    // calculate forces and potential
    atom **PermutationMap = NULL;
    MinimiseConstrainedPotential Minimiser(atoms);
    //double ConstrainedPotentialEnergy =
    Minimiser(PermutationMap, configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
    Minimiser.EvaluateConstrainedForces(&Force);
    delete[](PermutationMap);
  }

  // and perform Verlet integration for each atom with position, velocity and force vector
  // check size of vectors
  BOOST_FOREACH(atom *_atom, atoms) {
    _atom->VelocityVerletUpdate(_atom->getNr(), MDSteps+1, &configuration, &Force, (const size_t) 0);
  }

  // correct velocities (rather momenta) so that center of mass remains motionless
  Velocity = atoms.totalMomentumAtStep(MDSteps+1);
  IonMass = atoms.totalMass();

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

  atoms.addVelocityAtStep(-1*Velocity,MDSteps+1);
  ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1);
  Berendsen berendsen = Berendsen();
  berendsen.addToContainer(World::getInstance().getThermostats());
  double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
  MDSteps++;

  // exit
  return true;
};