/*
 * 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 <list>
#include <vector>

#include <gsl/gsl_randist.h>

#include "tesselation.hpp"

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

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

#define BondList list<bond *>

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

class AtomInfo {
public:
  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

  AtomInfo();
  ~AtomInfo();

private:
};

class TrajectoryParticleInfo {
public:
  struct
  {
    vector<Vector> R;  //!< position vector
    vector<Vector> U;  //!< velocity vector
    vector<Vector> F;  //!< last force vector
  } Trajectory;
  int FixedIon;   //!< config variable that states whether forces act on the ion or not

  TrajectoryParticleInfo();
  ~TrajectoryParticleInfo();

private:
};

class TrajectoryParticle :  public TrajectoryParticleInfo, public virtual AtomInfo, public virtual ParticleInfo {
public:

  TrajectoryParticle();
  virtual ~TrajectoryParticle();

  // 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);

  // 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 );

  // 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);

private:

};

class GraphNode;

class GraphNodeInfo {
public:

  int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
  int *ComponentNr;//!< belongs to this non-separable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
  int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect non-separable 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()
  atom *Ancestor; //!< "Father" in Depth-First-Search

  GraphNodeInfo();
  ~GraphNodeInfo();
private:

};


class GraphNode : public GraphNodeInfo, public virtual ParticleInfo {
public:

  GraphNode();
  virtual ~GraphNode();

  void OutputGraphInfo(ofstream *out) const;
  void OutputComponentNumber(ofstream *out) const;

private:

};

class BondedParticleInfo {
public:
  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
  BondList ListOfBonds; //!< list of all bonds

  BondedParticleInfo();
  ~BondedParticleInfo();

private:

};

class BondedParticle : public BondedParticleInfo, public virtual ParticleInfo, public virtual AtomInfo {
public:
  BondedParticle();
  virtual ~BondedParticle();

  bool RegisterBond(bond *Binder);
  bool UnregisterBond(bond *Binder);
  void UnregisterAllBond();
  int CountBonds() const;
  int CorrectBondDegree(ofstream *out);
  bool OutputBondOfAtom(ofstream *out) const;
  void OutputAdjacency(ofstream *AdjacencyFile) const;
  void OutputOrder(ofstream *file);

private:

};

/** Single atom.
 * Class incorporates position, type
 */
class atom : public TesselPoint, public TrajectoryParticle, public GraphNode, public BondedParticle, public virtual ParticleInfo, public virtual AtomInfo {
  public:
    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
    int *sort;      //!< sort criteria

  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;
  void OutputMPQCLine(ofstream *out, Vector *center, int *AtomNo) const;

  void InitComponentNr();

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

  double DistanceToVector(Vector &origin);
  double DistanceSquaredToVector(Vector &origin);
  bool IsInParallelepiped(Vector offset, double *parallelepiped);

  private:
};

#endif /* ATOM_HPP_ */