/*
 * atom.hpp
 *
 *  Created on: Aug 3, 2009
 *      Author: heber
 */

#ifndef ATOM_HPP_
#define ATOM_HPP_

using namespace std;

/*********************************************** includes ***********************************/

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

#include <iostream>
#include <vector>

#include <gsl/gsl_randist.h>

#include "tesselation.hpp"

/****************************************** forward declarations *****************************/

class bond;
class config;
class element;
class ForceMatrix;
class Vector;

/********************************************** declarations *******************************/

/** Single atom.
 * Class incorporates position, type
 */
class atom : public TesselPoint {
  public:
    struct
    {
      vector<Vector> R;  //!< position vector
      vector<Vector> U;  //!< velocity vector
      vector<Vector> F;  //!< last force vector
    } Trajectory;

    Vector x;       //!< coordinate vector of atom, giving last position within cell
    Vector v;       //!< velocity vector of atom, giving last velocity within cell
    Vector F;       //!< Force vector of atom, giving last force within cell
    element *type;  //!< pointing to element
    atom *previous; //!< previous atom in molecule list
    atom *next;     //!< next atom in molecule list
    atom *father;   //!< In many-body bond order fragmentations points to originating atom
    atom *Ancestor; //!< "Father" in Depth-First-Search
    //char *Name;      //!< unique name used during many-body bond-order fragmentation, comes from TesselPoint
    int FixedIon;   //!< config variable that states whether forces act on the ion or not
    int *sort;      //!< sort criteria
    //int nr;         //!< continuous, unique number, comes from TesselPoint
    int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
    int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
    int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
    bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
    bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not

  atom();
  atom(class atom *pointer);
  virtual ~atom();

  bool Output(ofstream *out, int ElementNo, int AtomNo, const char *comment = NULL) const;
  bool Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment = NULL);
  bool OutputXYZLine(ofstream *out) const;
  bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
  bool OutputTrajectoryXYZ(ofstream *out, int step) const;
  bool OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const;

  void EqualsFather ( atom *ptr, atom **res );
  void CorrectFather();
  atom *GetTrueFather();
  bool Compare(const atom &ptr);

  // trajectory stuff
  void ResizeTrajectory(int MaxSteps);
  void CopyStepOnStep(int dest, int src);
  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );

  double DistanceToVector(Vector &origin);
  double DistanceSquaredToVector(Vector &origin);

  bool IsInParallelepiped(Vector offset, double *parallelepiped);

  // bond order stuff
  void OutputOrder(ofstream *file);

  // constraint potential and dynamics stuff
  void AddKineticToTemperature(double *temperature, int step) const;
  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
  void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);

  // thermostats
  void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
  void Thermostat_Gaussian_init(int Step, double *G, double *E);
  void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
  void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
  void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
  void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
  void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);


  ostream & operator << (ostream &ost);

  private:
};

ostream & operator << (ostream &ost, const atom &a);

#endif /* ATOM_HPP_ */