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

/*
 * atom_atominfo.cpp
 *
 *  Created on: Oct 19, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "config.hpp"
#include "element.hpp"
#include "parser.hpp"
#include "periodentafel.hpp"
#include "World.hpp"
#include "atom_atominfo.hpp"

/** Constructor of class AtomInfo.
 */
AtomInfo::AtomInfo() :
  AtomicElement(NULL),
  FixedIon(false)
{
  AtomicPosition.reserve(1);
  AtomicPosition.push_back(zeroVec);
  AtomicVelocity.reserve(1);
  AtomicVelocity.push_back(zeroVec);
  AtomicForce.reserve(1);
  AtomicForce.push_back(zeroVec);
};

/** Copy constructor of class AtomInfo.
 */
AtomInfo::AtomInfo(const AtomInfo &_atom) :
    AtomicPosition(_atom.AtomicPosition),
    AtomicElement(_atom.AtomicElement),
    FixedIon(false)
{
  AtomicVelocity.reserve(1);
  AtomicVelocity.push_back(zeroVec);
  AtomicForce.reserve(1);
  AtomicForce.push_back(zeroVec);
};

AtomInfo::AtomInfo(const VectorInterface &_v) :
    AtomicElement(NULL),
    FixedIon(false)
{
  AtomicPosition[0] = _v.getPosition();
  AtomicVelocity.reserve(1);
  AtomicVelocity.push_back(zeroVec);
  AtomicForce.reserve(1);
  AtomicForce.push_back(zeroVec);
};

/** Destructor of class AtomInfo.
 */
AtomInfo::~AtomInfo()
{
};

const element *AtomInfo::getType() const
{
  return AtomicElement;
}

const double& AtomInfo::operator[](size_t i) const
{
  return AtomicPosition[0][i];
}

const double& AtomInfo::at(size_t i) const
{
  return AtomicPosition[0].at(i);
}

void AtomInfo::set(size_t i, const double value)
{
  AtomicPosition[0].at(i) = value;
}

const Vector& AtomInfo::getPosition() const
{
  return AtomicPosition[0];
}

const Vector& AtomInfo::getPosition(const int _step) const
{
  ASSERT(_step < AtomicPosition.size(),
      "AtomInfo::getPosition() - Access out of range.");
  return AtomicPosition[_step];
}

void AtomInfo::setType(const element* _type) {
  AtomicElement = _type;
}

void AtomInfo::setType(const int Z) {
  const element *elem = World::getInstance().getPeriode()->FindElement(Z);
  setType(elem);
}

Vector& AtomInfo::getAtomicVelocity()
{
  return AtomicVelocity[0];
}

Vector& AtomInfo::getAtomicVelocity(const int _step)
{
  ASSERT(_step < AtomicVelocity.size(),
      "AtomInfo::getAtomicVelocity() - Access out of range.");
  return AtomicVelocity[_step];
}

const Vector& AtomInfo::getAtomicVelocity() const
{
  return AtomicVelocity[0];
}

const Vector& AtomInfo::getAtomicVelocity(const int _step) const
{
  ASSERT(_step < AtomicVelocity.size(),
      "AtomInfo::getAtomicVelocity() - Access out of range.");
  return AtomicVelocity[_step];
}

void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
{
  AtomicVelocity[0] = _newvelocity;
}

void AtomInfo::setAtomicVelocity(const int _step, const Vector &_newvelocity)
{
  ASSERT(_step <= AtomicVelocity.size(),
      "AtomInfo::setAtomicVelocity() - Access out of range.");
  if(_step < (int)AtomicVelocity.size()) {
    AtomicVelocity[_step] = _newvelocity;
  } else if (_step == (int)AtomicVelocity.size()) {
    AtomicVelocity.push_back(_newvelocity);
  }
}

const Vector& AtomInfo::getAtomicForce() const
{
  return AtomicForce[0];
}

const Vector& AtomInfo::getAtomicForce(const int _step) const
{
  ASSERT(_step < AtomicForce.size(),
      "AtomInfo::getAtomicForce() - Access out of range.");
  return AtomicForce[_step];
}

void AtomInfo::setAtomicForce(const Vector &_newforce)
{
  AtomicForce[0] = _newforce;
}

void AtomInfo::setAtomicForce(const int _step, const Vector &_newforce)
{
  ASSERT(_step <= AtomicForce.size(),
      "AtomInfo::setAtomicForce() - Access out of range.");
  if(_step < (int)AtomicForce.size()) {
    AtomicForce[_step] = _newforce;
  } else if (_step == (int)AtomicForce.size()) {
    AtomicForce.push_back(_newforce);
  }
}

bool AtomInfo::getFixedIon() const
{
  return FixedIon;
}

void AtomInfo::setFixedIon(const bool _fixedion)
{
  FixedIon = _fixedion;
}

void AtomInfo::setPosition(const Vector& _vector)
{
  AtomicPosition[0] = _vector;
  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
}

void AtomInfo::setPosition(int _step, const Vector& _vector)
{
  ASSERT(_step <= AtomicPosition.size(),
      "AtomInfo::setPosition() - Access out of range.");
  if(_step < (int)AtomicPosition.size()) {
    AtomicPosition[_step] = _vector;
  } else if (_step == (int)AtomicPosition.size()) {
    AtomicPosition.push_back(_vector);
  }
  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
}

const VectorInterface& AtomInfo::operator+=(const Vector& b)
{
  AtomicPosition[0] += b;
  return *this;
}

const VectorInterface& AtomInfo::operator-=(const Vector& b)
{
  AtomicPosition[0] -= b;
  return *this;
}

Vector const AtomInfo::operator+(const Vector& b) const
{
  Vector a(AtomicPosition[0]);
  a += b;
  return a;
}

Vector const AtomInfo::operator-(const Vector& b) const
{
  Vector a(AtomicPosition[0]);
  a -= b;
  return a;
}

double AtomInfo::distance(const Vector &point) const
{
  return AtomicPosition[0].distance(point);
}

double AtomInfo::DistanceSquared(const Vector &y) const
{
  return AtomicPosition[0].DistanceSquared(y);
}

double AtomInfo::distance(const VectorInterface &_atom) const
{
  return _atom.distance(AtomicPosition[0]);
}

double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
{
  return _atom.DistanceSquared(AtomicPosition[0]);
}

VectorInterface &AtomInfo::operator=(const Vector& _vector)
{
  AtomicPosition[0] = _vector;
  return *this;
}

void AtomInfo::ScaleAll(const double *factor)
{
  AtomicPosition[0].ScaleAll(factor);
}

void AtomInfo::ScaleAll(const Vector &factor)
{
  AtomicPosition[0].ScaleAll(factor);
}

void AtomInfo::Scale(const double factor)
{
  AtomicPosition[0].Scale(factor);
}

void AtomInfo::Zero()
{
  AtomicPosition[0].Zero();
}

void AtomInfo::One(const double one)
{
  AtomicPosition[0].One(one);
}

void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
{
  AtomicPosition[0].LinearCombinationOfVectors(x1,x2,x3,factors);
}

/**
 *  returns the kinetic energy of this atom at a given time step
 */
double AtomInfo::getKineticEnergy(unsigned int step) const{
  return getMass() * AtomicVelocity[step].NormSquared();
}

Vector AtomInfo::getMomentum(unsigned int step) const{
  return getMass()*AtomicVelocity[step];
}

/** Extends the trajectory STL vector to the new size.
 * Does nothing if \a MaxSteps is smaller than current size.
 * \param MaxSteps
 */
void AtomInfo::ResizeTrajectory(size_t MaxSteps)
{
  if (AtomicPosition.size() <= (unsigned int)(MaxSteps)) {
    DoLog(0) && (Log() << Verbose(0) << "Increasing size for trajectory array from " << AtomicPosition.size() << " to " << (MaxSteps+1) << "." << endl);
    AtomicPosition.resize(MaxSteps+1, zeroVec);
    AtomicVelocity.resize(MaxSteps+1, zeroVec);
    AtomicForce.resize(MaxSteps+1, zeroVec);
  }
};

size_t AtomInfo::getTrajectorySize() const
{
  return AtomicPosition.size();
}

double AtomInfo::getMass() const{
  return AtomicElement->getMass();
}

/** Copies a given trajectory step \a src onto another \a dest
 * \param dest index of destination step
 * \param src index of source step
 */
void AtomInfo::CopyStepOnStep(int dest, int src)
{
  if (dest == src)  // self assignment check
    return;

  ASSERT(dest < AtomicPosition.size(),
      "AtomInfo::CopyStepOnStep() - destination outside of current trajectory array.");
  ASSERT(src < AtomicPosition.size(),
      "AtomInfo::CopyStepOnStep() - source outside of current trajectory array.");
  for (int n=NDIM;n--;) {
    AtomicPosition.at(dest)[n] = AtomicPosition.at(src)[n];
    AtomicVelocity.at(dest)[n] = AtomicVelocity.at(src)[n];
    AtomicForce.at(dest)[n] = AtomicForce.at(src)[n];
  }
};

/** Performs a velocity verlet update of the trajectory.
 * Parameters are according to those in configuration class.
 * \param NextStep index of sequential step to set
 * \param *configuration pointer to configuration with parameters
 * \param *Force matrix with forces
 */
void AtomInfo::VelocityVerletUpdate(int nr, int NextStep, config *configuration, ForceMatrix *Force, const size_t offset)
{
  //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
  if ((int)AtomicPosition.size() <= NextStep)
    ResizeTrajectory(NextStep+10);
  for (int d=0; d<NDIM; d++) {
    AtomicForce.at(NextStep)[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
    AtomicPosition.at(NextStep)[d] = AtomicPosition.at(NextStep-1)[d];
    AtomicPosition.at(NextStep)[d] += configuration->Deltat*(AtomicVelocity.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
    AtomicPosition.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(AtomicForce.at(NextStep)[d]/getMass());     // F = m * a and s =
  }
  // Update U
  for (int d=0; d<NDIM; d++) {
    AtomicVelocity.at(NextStep)[d] = AtomicVelocity.at(NextStep-1)[d];
    AtomicVelocity.at(NextStep)[d] += configuration->Deltat * (AtomicForce.at(NextStep)[d]+AtomicForce.at(NextStep-1)[d]/getMass()); // v = F/m * t
  }
  // Update R (and F)
//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
//      for (int d=0;d<NDIM;d++)
//        out << AtomicPosition.at(NextStep).x[d] << " ";          // next step
//      out << ")\t(";
//      for (int d=0;d<NDIM;d++)
//        Log() << Verbose(0) << AtomicVelocity.at(NextStep).x[d] << " ";          // next step
//      out << ")" << endl;
};

const AtomInfo& operator*=(AtomInfo& a, const double m)
{
  a.Scale(m);
  return a;
}

AtomInfo const operator*(const AtomInfo& a, const double m)
{
  AtomInfo copy(a);
  copy *= m;
  return copy;
}

AtomInfo const operator*(const double m, const AtomInfo& a)
{
  AtomInfo copy(a);
  copy *= m;
  return copy;
}

std::ostream & AtomInfo::operator << (std::ostream &ost) const
{
  return (ost << getPosition());
}

std::ostream & operator << (std::ostream &ost, const AtomInfo &a)
{
  ost << a;
  return ost;
}