/** \file molecules.hpp
 *
 * Class definitions of atom and molecule, element and periodentafel 
 */

#ifndef MOLECULES_HPP_
#define MOLECULES_HPP_

using namespace std;

// GSL headers
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_heapsort.h>

// STL headers
#include <map>
#include <set>
#include <deque>

#include "helpers.hpp"

class atom;
class AtomStackClass;
class bond;
class config;
class element;
class molecule;
class MoleculeListClass;
class periodentafel;
class vector;
class Verbose;

/******************************** Some definitions for easier reading **********************************/

#define KeyStack deque<int>
#define KeySet set<int>
#define Graph map<KeySet, pair<int, double>, KeyCompare >
#define GraphPair pair<KeySet, pair<int, double> >
#define KeySetTestPair pair<KeySet::iterator, bool>
#define GraphTestPair pair<Graph::iterator, bool>

struct KeyCompare
{
  bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
};
//bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)
inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment); // Insert a KeySet into a Graph
inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);  // Insert all KeySet's in a Graph into another Graph 

/******************************** Some templates for list management ***********************************/

/** Adds linking of an item to a list.
 * \param *walker
 * \return true - adding succeeded, false - error in list
 */
template <typename X> void link(X *walker, X *end)
{
  X *vorher = end->previous;
  if (vorher != NULL)
    vorher->next = walker;
  end->previous = walker;
  walker->previous = vorher;
  walker->next = end;
};

/** Removes linking of an item in a list.
 * \param *walker
 * \return true - removing succeeded, false - given item not found in list
 */
template <typename X> void unlink(X *walker)
{
  if (walker->next != NULL) 
    walker->next->previous = walker->previous;
  if (walker->previous != NULL)
    walker->previous->next = walker->next;
};

/** Adds new item before an item \a *end in a list.
 * \param *pointer   item to be added
 * \param *end  end of list
 * \return true - addition succeeded, false - unable to add item to list
 */
template <typename X>  bool add(X *pointer, X *end)
{
  if (end != NULL) {
    link(pointer, end);
  } else {
    pointer->previous = NULL;
    pointer->next = NULL; 
  }
  return true;
};

/** Finds item in list
 * \param *suche  search criteria
 * \param *start  begin of list
 * \param *end  end of list
 * \return X - if found, NULL - if not found
 */
template <typename X, typename Y> X * find(Y *suche, X *start, X *end)
{
  X *walker = start;
  while (walker->next != end) { // go through list
    walker = walker->next; // step onward beforehand
    if (*walker->sort == *suche) return (walker);
  }
  return NULL;
};

/** Removes an item from the list without check.
 * \param *walker item to be removed
 * \return true - removing succeeded, false - given item not found in list
 */
template <typename X> void removewithoutcheck(X *walker)
{
  if (walker != NULL) {
    unlink(walker);
    delete(walker);
    walker = NULL;
  }
};

/** Removes an item from the list, checks if exists.
 * Checks beforehand if atom is really within molecule list.
 * \param *pointer   item to be removed
 * \param *start  begin of list
 * \param *end  end of list
 * \return true - removing succeeded, false - given item not found in list
 */
template <typename X> bool remove(X *pointer, X *start, X *end)
{
  X *walker = find (pointer->sort, start, end);
/*  while (walker->next != pointer) { // search through list
    walker = walker->next;
    if (walker == end) return false;  // item not found in list
  }*/
  // atom found, now unlink
  if (walker != NULL)
    removewithoutcheck(walker);
  else
    return false;
  return true;
};

/** Cleans the whole list.
 * \param *start begin of list
 * \param *end end of list
 * \return true - list was cleaned successfully, false - error in list structure
 */
template <typename X> bool cleanup(X *start, X *end)
{
  X *pointer = start->next;
  X *walker;
  while (pointer != end) { // go through list
    walker = pointer; // mark current
    pointer = pointer->next; // step onward beforehand
    // remove walker
    unlink(walker);
    delete(walker);
    walker = NULL;
  }
  return true;
};

/** Returns the first marker in a chain list.
 * \param *me one arbitrary item in chain list
 * \return poiner to first marker
 */
template <typename X> X *GetFirst(X *me)
{
  X *Binder = me;
  while(Binder->previous != NULL)
    Binder = Binder->previous;
  return Binder;
}; 

/** Returns the last marker in a chain list.
 * \param *me one arbitrary item in chain list
 * \return poiner to last marker
 */
template <typename X> X *GetLast(X *me)
{
  X *Binder = me;
  while(Binder->next != NULL)
    Binder = Binder->next;
  return Binder;
}; 

/** Frees a two-dimensional array.
 * \param *ptr pointer to array
 * \param dim first dim of array
 */
template <typename X> void Free2DArray(X **ptr, int dim)
{
	int i;
 	if (ptr != NULL) {
 		for(i=0;i<dim;i++)
 			if (ptr[i] != NULL)
 				free(ptr[i]);
 		free(ptr);
 	}
};

int CompareDoubles (const void * a, const void * b);

/************************************* Class definitions ****************************************/

/** Chemical element.
 * Class incorporates data for a certain chemical element to be referenced from atom class.
 */
class element {
  public:
    double mass;    //!< mass in g/mol
    double CovalentRadius;  //!< covalent radius
    double VanDerWaalsRadius;  //!< can-der-Waals radius
    int Z;          //!< atomic number
    char name[64];  //!< atom name, i.e. "Hydrogren"
    char symbol[3]; //!< short form of the atom, i.e. "H"
    char period[8];    //!< period: n quantum number
    char group[8];    //!< group: l quantum number
    char block[8];    //!< block: l quantum number
    element *previous;  //!< previous item in list
    element *next;  //!< next element in list
    int *sort;      //!< sorc criteria
    int No;         //!< number of element set on periodentafel::Output()
    double Valence;		//!< number of valence electrons for this element
    int NoValenceOrbitals;  //!< number of valence orbitals, used for determining bond degree in molecule::CreateConnectmatrix()
    double HBondDistance[NDIM];	//!< distance in Angstrom of this element to hydrogen  (for single, double and triple bonds)
    double HBondAngle[NDIM];     //!< typical angle for one, two, three bonded hydrogen (in degrees)

  element();
  ~element();

  //> print element entries to screen
  bool Output(ofstream *out) const;
  bool Checkout(ofstream *out, const int No, const int NoOfAtoms) const;
  
  private:
};

/** Periodentafel is a list of all elements sorted by their atomic number.
 */
class periodentafel {
  public:
    element *start; //!< start of element list
    element *end;   //!< end of element list
    char header1[255]; //!< store first header line
    char header2[255]; //!< store second header line
  
  periodentafel();
  ~periodentafel(); 
  
  bool AddElement(element *pointer);
  bool RemoveElement(element *pointer);
  bool CleanupPeriodtable();
  element * FindElement(int Z);
  element * FindElement(char *shorthand) const;
  element * AskElement(); 
  bool Output(ofstream *output) const;
  bool Checkout(ofstream *output, const int *checkliste) const;
  bool LoadPeriodentafel(char *filename = NULL);
  bool StorePeriodentafel(char *filename = NULL) const;
  
  private:
};

// some algebraic matrix stuff
#define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3])  //!< hard-coded determinant of a 3x3 matrix
#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                      //!< hard-coded determinant of a 2x2 matrix

/** Single vector.
 * basically, just a x[3] but with helpful functions
 */
class vector {   
  public:
    double x[NDIM];

  vector();
  ~vector();

  double Distance(const vector *y) const;
  double PeriodicDistance(const vector *y, const double *cell_size) const;
  double ScalarProduct(const vector *y) const;
  double Projection(const vector *y) const;
  double Norm() const ;
  double Angle(vector *y) const;

  void AddVector(const vector *y);
  void SubtractVector(const vector *y);
  void CopyVector(const vector *y);
  void RotateVector(const vector *y, const double alpha);
  void Zero(); 
  void Normalize();
  void Translate(const vector *x);
  void Mirror(const vector *x);
  void Scale(double **factor);
  void Scale(double *factor);
  void Scale(double factor);
  void MatrixMultiplication(double *M);
  void InverseMatrixMultiplication(double *M);
  void KeepPeriodic(ofstream *out, double *matrix);
  void LinearCombinationOfVectors(const vector *x1, const vector *x2, const vector *x3, double *factors);
  
  bool GetOneNormalVector(const vector *x1);
  bool MakeNormalVector(const vector *y1);
  bool MakeNormalVector(const vector *y1, const vector *y2);
  bool MakeNormalVector(const vector *x1, const vector *x2, const vector *x3);
  bool SolveSystem(vector *x1, vector *x2, vector *y, double alpha, double beta, double c);
	bool LSQdistance(vector **vectors, int dim);

  void AskPosition(double *cell_size, bool check);
  bool Output(ofstream *out) const;
};

ofstream& operator<<(ofstream& ost, vector& m);

/** Parameter structure for least square minimsation.
 */
struct LSQ_params {
  vector **vectors;
  int num;
};

double LSQ(const gsl_vector * x, void * params);

/** Parameter structure for least square minimsation.
 */
struct lsq_params {
  gsl_vector *x;
  const molecule *mol;
  element *type;
};



/** Single atom.
 * Class incoporates position, type
 */
class atom {
  public:
    vector x;       //!< coordinate array of atom, giving position within cell
    vector v;       //!< velocity array of atom
    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
    int FixedIon;   //!< config variable that states whether forces act on the ion or not
    int *sort;      //!< sort criteria
    int nr;         //!< continuous, unique number
    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, given in DepthFirstSearchAnalysis()
  
  atom();
  ~atom();
  
  bool Output(int ElementNo, int AtomNo, ofstream *out) const;
  bool OutputXYZLine(ofstream *out) const;
  atom *GetTrueFather();
  bool Compare(atom &ptr);
  
  private:
};

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

/* Stack of Atoms.
 * Is used during DepthFirstSearchAnalysis() to detect nonseparable components.
 */
class AtomStackClass {
  public:
  AtomStackClass(int dimension);
  ~AtomStackClass();
  
  bool Push(atom *object);
  atom *PopFirst();
  atom *PopLast();
  bool AtomStackClass::RemoveItem(atom *ptr);
  void ClearStack();
  bool IsEmpty();
  bool IsFull();
	int ItemCount();
  void Output(ofstream *out) const;
  void TestImplementation(ofstream *out, atom *test);
  
  private:
    atom **StackList;   //!< the list containing the atom pointers
    int EntryCount;     //!< number of entries in the stack
    int CurrentLastEntry;   //!< Current last entry (newest item on stack)
    int CurrentFirstEntry;   //!< Current first entry (oldest item on stack)
    int NextFreeField;       //!< Current index of next free field
};


/** Bonds between atoms.
 * Class incorporates bonds between atoms in a molecule,
 * used to derive tge fragments in many-body bond order
 * calculations.
 */
class bond {
  public:
  	atom *leftatom;		//!< first bond partner
  	atom *rightatom;	//!< second bond partner
    bond *previous; //!< previous atom in molecule list
    bond *next;     //!< next atom in molecule list
  	int HydrogenBond;	//!< Number of hydrogen atoms in the bond
  	int BondDegree;		//!< single, double, triple, ... bond
    int nr;           //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
    bool Cyclic;      //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
    enum EdgeType Type;//!< whether this is a tree or back edge 
  	
  atom * GetOtherAtom(atom *Atom) const;
  bond * GetFirstBond();
  bond * GetLastBond();
  
  bool MarkUsed(enum Shading color);
  enum Shading IsUsed();
  void ResetUsed();
  bool Contains(const atom *ptr);
  bool Contains(const int nr);
  
  bond();
  bond(atom *left, atom *right);
  bond(atom *left, atom *right, int degree);
  bond(atom *left, atom *right, int degree, int number);
  ~bond();
    
  private: 
    enum Shading Used;        //!< marker in depth-first search, DepthFirstSearchAnalysis()
};

ostream & operator << (ostream &ost, bond &b);

class MoleculeLeafClass;

/** The complete molecule.
 * Class incorporates number of types
 */
class molecule {
  public:
    double cell_size[6];//!< cell size
    periodentafel *elemente; //!< periodic table with each element
    atom *start;        //!< start of atom list
    atom *end;          //!< end of atom list
    bond *first;        //!< start of bond list
    bond *last;         //!< end of bond list
    bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
    int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has
    int AtomCount;					//!< number of atoms, brought up-to-date by CountAtoms()
    int BondCount;					//!< number of atoms, brought up-to-date by CountBonds()
    int ElementCount;       //!< how many unique elements are therein
    int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
    int NoNonHydrogen;	//!< number of non-hydrogen atoms in molecule
    int NoNonBonds;     //!< number of non-hydrogen bonds in molecule
    int NoCyclicBonds;  //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
    double BondDistance;  //!< typical bond distance used in CreateAdjacencyList() and furtheron
  
  molecule(periodentafel *teil);
  ~molecule();
  
  /// remove atoms from molecule.
  bool AddAtom(atom *pointer);
  bool RemoveAtom(atom *pointer);
  bool CleanupMolecule();

  /// Add/remove atoms to/from molecule.
  atom * AddCopyAtom(atom *pointer);
  bool AddXYZFile(string filename);
  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);  
  bond * AddBond(atom *first, atom *second, int degree);
  bool RemoveBond(bond *pointer);
  bool RemoveBonds(atom *BondPartner);
    
  /// Find atoms.
  atom * FindAtom(int Nr) const; 
  atom * AskAtom(char *text);

  /// Count and change present atoms' coordination.
  void CountAtoms(ofstream *out);
  void CountElements();
  void CalculateOrbitals(class config &configuration);
  void CenterEdge(ofstream *out, vector *max); 
  void CenterOrigin(ofstream *out, vector *max); 
  void CenterGravity(ofstream *out, vector *max); 
  void Translate(const vector *x);
  void Mirror(const vector *x);
  void Align(vector *n);
  void Scale(double **factor);
  void DetermineCenterOfGravity(vector &CenterOfGravity);
  void SetBoxDimension(vector *dim);
  double * ReturnFullMatrixforSymmetric(double *cell_size);
  void ScanForPeriodicCorrection(ofstream *out);

  bool CheckBounds(const vector *x) const;
  void GetAlignVector(struct lsq_params * par) const;

  /// Initialising routines in fragmentation  
  void CreateAdjacencyList(ofstream *out, double bonddistance);
  void CreateListOfBondsPerAtom(ofstream *out);
  
  // Graph analysis
  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int &MinimumRingSize);
  void CyclicStructureAnalysis(ofstream *out, int &MinimumRingSize);
  bond * FindNextUnused(atom *vertex);
  void SetNextComponentNumber(atom *vertex, int nr);
  void InitComponentNumbers();
  void OutputComponentNumber(ofstream *out, atom *vertex);
  void ResetAllBondsToUnused();
  void ResetAllAtomNumbers();
  int CountCyclicBonds(ofstream *out);
  char * GetColor(enum Shading color);

  molecule *CopyMolecule();
  
  /// Fragment molecule by two different approaches:
  void FragmentMolecule(ofstream *out, int BottomUpOrder, int TopDownOrder, enum BondOrderScheme Scheme, config *configuration, enum CutCyclicBond CutCyclic);
  bool StoreAdjacencyToFile(ofstream *out, char *path);
  bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
  bool ParseKeySetFile(ofstream *out, char *filename, atom **ListOfAtoms, MoleculeListClass *&FragmentList, bool IsAngstroem);
  bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
  MoleculeListClass * GetAtomicFragments(ofstream *out, int NumberOfTopAtoms, bool IsAngstroem, double factor, enum CutCyclicBond CutCyclic);
  void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem, enum CutCyclicBond CutCyclic);
  /// -# BottomUp  
  MoleculeListClass * FragmentBottomUp(ofstream *out, int BondOrder, config *configuration, enum CutCyclicBond CutCyclic);
  MoleculeListClass * GetEachBondFragmentOfOrder(ofstream *out, int BondOrder, bool IsAngstroem, enum CutCyclicBond CutCyclic);
  /// -# TopDown
  void FragmentMoleculeByBond(ofstream *out, bond *Bond, molecule **LeftFragment, molecule **RightFragment, bool IsAngstroem, enum CutCyclicBond CutCyclic);
  /// -# BOSSANOVA
  Graph * FragmentBOSSANOVA(ofstream *out, int ANOVAOrder, config *configuration);
	int CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, Graph *ListOfGraph, KeySet Fragment, double TEFactor, config *configuration);
  bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
  molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
  void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
  int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
  int GuesstimateFragmentCount(ofstream *out, int order);
  	  
  // Recognize doubly appearing molecules in a list of them   
  int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
  int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);
	
  // Output routines.
  bool Output(ofstream *out);
  bool OutputXYZ(ofstream *out) const;
  bool Checkout(ofstream *out) const;

  private:
  int last_atom;      //!< number given to last atom
};

/** A list of \a molecule classes.
 */
class MoleculeListClass {
  public:
    molecule **ListOfMolecules;   //!< pointer list of fragment molecules to check for equality
    int NumberOfMolecules;        //!< Number of entries in \a **FragmentList and of to be returned one.
    int NumberOfTopAtoms;         //!< Number of atoms in the molecule from which all fragments originate
    double *TEList;               //!< List of factors when summing over total energies of all fragment
    
  MoleculeListClass();
  MoleculeListClass(int Num, int NumAtoms);
  ~MoleculeListClass();

  MoleculeListClass * FragmentTopDown(ofstream *out, int BondDegree, double bonddistance, config *configuration, enum CutCyclicBond Saturation);

  /// Reduced list to unique molecules.
  int * GetMappingToUniqueFragments(ofstream *out, double threshold, double *cell_size, double celldistance);
  int ReduceFragmentToUniqueOnes(ofstream *out, int *Map);
  void ReduceToUniqueList(ofstream *out, double *cell_size, double celldistance);

  /// Output configs.
  bool OutputConfigForListOfFragments(char *prefix, config *configuration, int *SortIndex);
  void Output(ofstream *out);
  
  private:
};


/** A leaf for a tree of \a molecule class
 * Wraps molecules in a tree structure
 */
class MoleculeLeafClass {
  public:
    molecule *Leaf;                   //!< molecule of this leaf 
    //MoleculeLeafClass *UpLeaf;        //!< Leaf one level up
    //MoleculeLeafClass *DownLeaf;      //!< First leaf one level down
    MoleculeLeafClass *previous;  //!< Previous leaf on this level
    MoleculeLeafClass *next;      //!< Next leaf on this level

  //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
  MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
  ~MoleculeLeafClass();

  bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
};

/** The config file.
 * The class contains all parameters that control a dft run also functions to load and save.
 */
class config {
  public:
    int PsiType;
    int MaxPsiDouble;
    int PsiMaxNoUp;
    int PsiMaxNoDown;
    int MaxMinStopStep;
    int InitMaxMinStopStep;
    int ProcPEGamma;
    int ProcPEPsi;
    char *configpath;
    char *configname;
    
    private:
    char *mainname;
    char *defaultpath;
    char *pseudopotpath;
    
    int DoOutVis;
    int DoOutMes;
    int DoOutNICS;
    int DoOutOrbitals;
    int DoOutCurrent;
    int DoFullCurrent;
    int DoPerturbation;
    int CommonWannier;
    double SawtoothStart;
    int VectorPlane;
    double VectorCut;
    int UseAddGramSch;
    int Seed;
    
    int MaxOuterStep;
    double Deltat;
    int OutVisStep;
    int OutSrcStep;
    double TargetTemp;
    int ScaleTempStep;
    int MaxPsiStep;
    double EpsWannier;
    
    int MaxMinStep;
    double RelEpsTotalEnergy;
    double RelEpsKineticEnergy;
    int MaxMinGapStopStep;
    int MaxInitMinStep;
    double InitRelEpsTotalEnergy;
    double InitRelEpsKineticEnergy;
    int InitMaxMinGapStopStep;
    
    //double BoxLength[NDIM*NDIM];
    
    double ECut;
    int MaxLevel;
    int RiemannTensor;
    int LevRFactor;
    int RiemannLevel;
    int Lev0Factor;
    int RTActualUse;
    int AddPsis;
    
    double RCut;
    int StructOpt;
    int IsAngstroem;
    int RelativeCoord;
    int MaxTypes;

  
  int ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical);
  
  public:
  config();
  ~config();

  int TestSyntax(char *filename, periodentafel *periode, molecule *mol);
  void Load(char *filename, periodentafel *periode, molecule *mol);
  void LoadOld(char *filename, periodentafel *periode, molecule *mol);
  void RetrieveConfigPathAndName(char *filename);
  bool Save(ofstream *file, periodentafel *periode, molecule *mol) const;
  void Edit(molecule *mol);
  bool GetIsAngstroem() const;
  char *GetDefaultPath() const;
  void config::SetDefaultPath(const char *path);
};

#endif /*MOLECULES_HPP_*/