source: src/molecules.hpp@ 6ac7ee

/** \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_eigen.h>
#include <gsl/gsl_heapsort.h>

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

#include "helpers.hpp"
#include "parser.hpp"
#include "periodentafel.hpp"
#include "stackclass.hpp"
#include "vector.hpp"

class atom;
class bond;
class config;
class molecule;
class MoleculeListClass;
class Verbose;

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

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

#define DistancePair pair < double, atom* >
#define DistanceMap multimap < double, atom* >
#define DistanceTestPair pair < DistanceMap::iterator, bool>

#define Boundaries map <double, DistancePair >
#define BoundariesPair pair<double, DistancePair >
#define BoundariesTestPair pair< Boundaries::iterator, bool>

#define PointMap map < int, class BoundaryPointSet * >
#define PointPair pair < int, class BoundaryPointSet * >
#define PointTestPair pair < PointMap::iterator, bool >

#define LineMap map < int, class BoundaryLineSet * >
#define LinePair pair < int, class BoundaryLineSet * >
#define LineTestPair pair < LineMap::iterator, bool >

#define TriangleMap map < int, class BoundaryTriangleSet * >
#define TrianglePair pair < int, class BoundaryTriangleSet * >
#define TriangleTestPair pair < TrianglePair::iterator, bool >

#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >

/******************************** Some small functions and/or structures **********************************/

struct KeyCompare
{
        bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
};

struct Trajectory
{
        vector<Vector> R;       //!< position vector
        vector<Vector> U;       //!< velocity vector
        vector<Vector> F;       //!< last force vector
        atom *ptr;                               //!< pointer to atom whose trajectory we contain
};

//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
int CompareDoubles (const void * a, const void * b);


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


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


/** 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, 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();

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

/** 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
                map<atom *, struct Trajectory> Trajectories; //!< contains old trajectory points
                int MDSteps;                            //!< The number of MD steps in Trajectories
                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(string text);

        /// Count and change present atoms' coordination.
        void CountAtoms(ofstream *out);
        void CountElements();
        void CalculateOrbitals(class config &configuration);
        bool CenterInBox(ofstream *out, Vector *BoxLengths);
        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 DetermineCenter(Vector &center);
        Vector * DetermineCenterOfGravity(ofstream *out);
        Vector * DetermineCenterOfAll(ofstream *out);
        void SetBoxDimension(Vector *dim);
        double * ReturnFullMatrixforSymmetric(double *cell_size);
        void ScanForPeriodicCorrection(ofstream *out);
        void PrincipalAxisSystem(ofstream *out, bool DoRotate);
        double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
        bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);

        bool CheckBounds(const Vector *x) const;
        void GetAlignvector(struct lsq_params * par) const;

        /// Initialising routines in fragmentation
        void CreateAdjacencyList2(ofstream *out, ifstream *output);
        void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
        void CreateListOfBondsPerAtom(ofstream *out);

        // Graph analysis
        MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack);
        void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
        bool PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack);
        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);
        bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
        string GetColor(enum Shading color);

        molecule *CopyMolecule();

        /// Fragment molecule by two different approaches:
        int FragmentMolecule(ofstream *out, int Order, config *configuration);
        bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
        bool StoreAdjacencyToFile(ofstream *out, char *path);
        bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
        bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
        bool StoreOrderAtSiteFile(ofstream *out, char *path);
        bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
        bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
        bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
        bool CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex);
        bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
        void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
        /// -# BOSSANOVA
        void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
        int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
        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 OutputTrajectories(ofstream *out);
        void OutputListOfBonds(ofstream *out) const;
        bool OutputXYZ(ofstream *out) const;
        bool OutputTrajectoriesXYZ(ofstream *out);
        bool Checkout(ofstream *out) const;
        bool OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output);

        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

        MoleculeListClass();
        MoleculeListClass(int Num, int NumAtoms);
        ~MoleculeListClass();

        /// Output configs.
        bool AddHydrogenCorrection(ofstream *out, char *path);
        bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
        bool OutputConfigForListOfFragments(ofstream *out, 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);
        bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
        bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
        bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
        bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList);
        void TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
        int Count() const;
};

/** 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;
                bool FastParsing;
                double Deltat;

                private:
                char *mainname;
                char *defaultpath;
                char *pseudopotpath;

                int DoOutVis;
                int DoOutMes;
                int DoOutNICS;
                int DoOutOrbitals;
                int DoOutCurrent;
                int DoFullCurrent;
                int DoPerturbation;
                int DoWannier;
                int CommonWannier;
                double SawtoothStart;
                int VectorPlane;
                double VectorCut;
                int UseAddGramSch;
                int Seed;

                int MaxOuterStep;
                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(string filename);
        bool Save(const char *filename, periodentafel *periode, molecule *mol) const;
        bool SaveMPQC(const char *filename, molecule *mol) const;
        void Edit(molecule *mol);
        bool GetIsAngstroem() const;
        char *GetDefaultPath() const;
        void SetDefaultPath(const char *path);
};

#endif /*MOLECULES_HPP_*/