Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/molecules.hpp

    redb93c ree1b16  
    1010
    1111// GSL headers
     12#include <gsl/gsl_multimin.h>
     13#include <gsl/gsl_vector.h>
     14#include <gsl/gsl_matrix.h>
    1215#include <gsl/gsl_eigen.h>
    1316#include <gsl/gsl_heapsort.h>
    14 #include <gsl/gsl_linalg.h>
    15 #include <gsl/gsl_matrix.h>
    16 #include <gsl/gsl_multimin.h>
    17 #include <gsl/gsl_vector.h>
    18 #include <gsl/gsl_randist.h>
    19 
    20 //// STL headers
     17
     18// STL headers
    2119#include <map>
    2220#include <set>
     
    2523#include <vector>
    2624
    27 #include "atom.hpp"
    28 #include "bond.hpp"
    29 #include "element.hpp"
    30 #include "linkedcell.hpp"
     25#include "helpers.hpp"
    3126#include "parser.hpp"
    3227#include "periodentafel.hpp"
    3328#include "stackclass.hpp"
    34 #include "tesselation.hpp"
    3529#include "vector.hpp"
    3630
     31class atom;
     32class bond;
     33class config;
    3734class molecule;
    38 class MoleculeLeafClass;
    3935class MoleculeListClass;
     36class Verbose;
    4037
    4138/******************************** Some definitions for easier reading **********************************/
     
    4946#define GraphTestPair pair<Graph::iterator, bool>
    5047
    51 #define MoleculeList list <molecule *>
    52 #define MoleculeListTest pair <MoleculeList::iterator, bool>
    53 
    5448#define DistancePair pair < double, atom* >
    5549#define DistanceMap multimap < double, atom* >
    5650#define DistanceTestPair pair < DistanceMap::iterator, bool>
    5751
    58 
    59 //#define LinkedAtoms list <atom *>
     52#define Boundaries map <double, DistancePair >
     53#define BoundariesPair pair<double, DistancePair >
     54#define BoundariesTestPair pair< Boundaries::iterator, bool>
     55
     56#define PointMap map < int, class BoundaryPointSet * >
     57#define PointPair pair < int, class BoundaryPointSet * >
     58#define PointTestPair pair < PointMap::iterator, bool >
     59#define CandidateList list <class CandidateForTesselation *>
     60
     61#define LineMap multimap < int, class BoundaryLineSet * >
     62#define LinePair pair < int, class BoundaryLineSet * >
     63#define LineTestPair pair < LineMap::iterator, bool >
     64
     65#define TriangleMap map < int, class BoundaryTriangleSet * >
     66#define TrianglePair pair < int, class BoundaryTriangleSet * >
     67#define TriangleTestPair pair < TrianglePair::iterator, bool >
     68
     69#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
     70#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
     71
     72#define MoleculeList list <molecule *>
     73#define MoleculeListTest pair <MoleculeList::iterator, bool>
     74
     75#define LinkedAtoms list <atom *>
    6076
    6177/******************************** Some small functions and/or structures **********************************/
     
    6379struct KeyCompare
    6480{
    65   bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
     81        bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
    6682};
    6783
    6884struct Trajectory
    6985{
    70   vector<Vector> R;  //!< position vector
    71   vector<Vector> U;  //!< velocity vector
    72   vector<Vector> F;  //!< last force vector
    73   atom *ptr;        //!< pointer to atom whose trajectory we contain
    74 };
    75 
    76 //bool operator < (KeySet SubgraphA, KeySet SubgraphB);  //note: this declaration is important, otherwise normal < is used (producing wrong order)
     86        vector<Vector> R;       //!< position vector
     87        vector<Vector> U;       //!< velocity vector
     88        vector<Vector> F;       //!< last force vector
     89        atom *ptr;                              //!< pointer to atom whose trajectory we contain
     90};
     91
     92//bool operator < (KeySet SubgraphA, KeySet SubgraphB);  //note: this declaration is important, otherwise normal < is used (producing wrong order)
    7793inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment); // Insert a KeySet into a Graph
    78 inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);  // Insert all KeySet's in a Graph into another Graph
     94inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);    // Insert all KeySet's in a Graph into another Graph
    7995int CompareDoubles (const void * a, const void * b);
    8096
     
    84100
    85101// some algebraic matrix stuff
    86 #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
    87 #define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                      //!< hard-coded determinant of a 2x2 matrix
     102#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
     103#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                                                                                        //!< hard-coded determinant of a 2x2 matrix
    88104
    89105
     
    91107 */
    92108struct LSQ_params {
    93   Vector **vectors;
    94   int num;
     109        Vector **vectors;
     110        int num;
    95111};
    96112
     
    100116 */
    101117struct lsq_params {
    102   gsl_vector *x;
    103   const molecule *mol;
    104   element *type;
    105 };
    106 
    107 
    108 #define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
    109 enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
    110 
     118        gsl_vector *x;
     119        const molecule *mol;
     120        element *type;
     121};
     122
     123/** Single atom.
     124 * Class incoporates position, type
     125 */
     126class atom {
     127        public:
     128                Vector x;                        //!< coordinate array of atom, giving position within cell
     129                Vector v;                        //!< velocity array of atom
     130                element *type;  //!< pointing to element
     131                atom *previous; //!< previous atom in molecule list
     132                atom *next;              //!< next atom in molecule list
     133                atom *father;    //!< In many-body bond order fragmentations points to originating atom
     134                atom *Ancestor; //!< "Father" in Depth-First-Search
     135                char *Name;                     //!< unique name used during many-body bond-order fragmentation
     136                int FixedIon;    //!< config variable that states whether forces act on the ion or not
     137                int *sort;                      //!< sort criteria
     138                int nr;                          //!< continuous, unique number
     139                int GraphNr;                    //!< unique number, given in DepthFirstSearchAnalysis()
     140                int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
     141                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.
     142                bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
     143                bool IsCyclic;                          //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
     144                unsigned char AdaptiveOrder;    //!< current present bond order at site (0 means "not set")
     145                bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
     146
     147        atom();
     148        ~atom();
     149
     150        bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
     151        bool OutputXYZLine(ofstream *out) const;
     152        atom *GetTrueFather();
     153        bool Compare(const atom &ptr);
     154
     155        private:
     156};
     157
     158ostream & operator << (ostream &ost, const atom &a);
     159
     160/** Bonds between atoms.
     161 * Class incorporates bonds between atoms in a molecule,
     162 * used to derive tge fragments in many-body bond order
     163 * calculations.
     164 */
     165class bond {
     166        public:
     167                atom *leftatom;         //!< first bond partner
     168                atom *rightatom;        //!< second bond partner
     169                bond *previous; //!< previous atom in molecule list
     170                bond *next;              //!< next atom in molecule list
     171                int HydrogenBond;       //!< Number of hydrogen atoms in the bond
     172                int BondDegree;         //!< single, double, triple, ... bond
     173                int nr;                                  //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
     174                bool Cyclic;                    //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
     175                enum EdgeType Type;//!< whether this is a tree or back edge
     176
     177        atom * GetOtherAtom(atom *Atom) const;
     178        bond * GetFirstBond();
     179        bond * GetLastBond();
     180
     181        bool MarkUsed(enum Shading color);
     182        enum Shading IsUsed();
     183        void ResetUsed();
     184        bool Contains(const atom *ptr);
     185        bool Contains(const int nr);
     186
     187        bond();
     188        bond(atom *left, atom *right);
     189        bond(atom *left, atom *right, int degree);
     190        bond(atom *left, atom *right, int degree, int number);
     191        ~bond();
     192
     193        private:
     194                enum Shading Used;                              //!< marker in depth-first search, DepthFirstSearchAnalysis()
     195};
     196
     197ostream & operator << (ostream &ost, bond &b);
     198
     199class MoleculeLeafClass;
    111200
    112201/** The complete molecule.
    113202 * Class incorporates number of types
    114203 */
    115 class molecule : public PointCloud {
    116   public:
    117     double cell_size[6];//!< cell size
    118     periodentafel *elemente; //!< periodic table with each element
    119     atom *start;        //!< start of atom list
    120     atom *end;          //!< end of atom list
    121     bond *first;        //!< start of bond list
    122     bond *last;         //!< end of bond list
    123     bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
    124     map<atom *, struct Trajectory> Trajectories; //!< contains old trajectory points
    125     int MDSteps;        //!< The number of MD steps in Trajectories
    126     int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has
    127     int AtomCount;          //!< number of atoms, brought up-to-date by CountAtoms()
    128     int BondCount;          //!< number of atoms, brought up-to-date by CountBonds()
    129     int ElementCount;       //!< how many unique elements are therein
    130     int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
    131     int NoNonHydrogen;  //!< number of non-hydrogen atoms in molecule
    132     int NoNonBonds;     //!< number of non-hydrogen bonds in molecule
    133     int NoCyclicBonds;  //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
    134     double BondDistance;  //!< typical bond distance used in CreateAdjacencyList() and furtheron
    135     bool ActiveFlag;    //!< in a MoleculeListClass used to discern active from inactive molecules
    136     Vector Center;      //!< Center of molecule in a global box
    137     char name[MAXSTRINGSIZE];         //!< arbitrary name
    138     int IndexNr;        //!< index of molecule in a MoleculeListClass
    139     class Tesselation *TesselStruct;
    140 
    141   molecule(periodentafel *teil);
    142   ~molecule();
    143 
    144   // re-definition of virtual functions from PointCloud
    145   Vector *GetCenter(ofstream *out);
    146   TesselPoint *GetPoint();
    147   TesselPoint *GetTerminalPoint();
    148   void GoToNext();
    149   void GoToPrevious();
    150   void GoToFirst();
    151   void GoToLast();
    152   bool IsEmpty();
    153   bool IsLast();
    154 
    155   /// remove atoms from molecule.
    156   bool AddAtom(atom *pointer);
    157   bool RemoveAtom(atom *pointer);
    158   bool UnlinkAtom(atom *pointer);
    159   bool CleanupMolecule();
    160 
    161   /// Add/remove atoms to/from molecule.
    162   atom * AddCopyAtom(atom *pointer);
    163   bool AddXYZFile(string filename);
    164   bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
    165   bond * AddBond(atom *first, atom *second, int degree);
    166   bool RemoveBond(bond *pointer);
    167   bool RemoveBonds(atom *BondPartner);
    168 
    169   /// Find atoms.
    170   atom * FindAtom(int Nr) const;
    171   atom * AskAtom(string text);
    172 
    173   /// Count and change present atoms' coordination.
    174   void CountAtoms(ofstream *out);
    175   void CountElements();
    176   void CalculateOrbitals(class config &configuration);
    177   bool CenterInBox(ofstream *out);
    178   void CenterEdge(ofstream *out, Vector *max);
    179   void CenterOrigin(ofstream *out);
    180   void CenterPeriodic(ofstream *out);
    181   void CenterAtVector(ofstream *out, Vector *newcenter);
    182   void Translate(const Vector *x);
    183   void TranslatePeriodically(const Vector *trans);
    184   void Mirror(const Vector *x);
    185   void Align(Vector *n);
    186   void Scale(double **factor);
    187   void DeterminePeriodicCenter(Vector &center);
    188   Vector * DetermineCenterOfGravity(ofstream *out);
    189   Vector * DetermineCenterOfAll(ofstream *out);
    190   void SetNameFromFilename(const char *filename);
    191   void SetBoxDimension(Vector *dim);
    192   double * ReturnFullMatrixforSymmetric(double *cell_size);
    193   void ScanForPeriodicCorrection(ofstream *out);
    194   bool VerletForceIntegration(ofstream *out, char *file, config &configuration);
    195   void Thermostats(config &configuration, double ActualTemp, int Thermostat);
    196   void PrincipalAxisSystem(ofstream *out, bool DoRotate);
    197   double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
    198   Vector* FindEmbeddingHole(ofstream *out, molecule *srcmol);
    199 
    200 
    201   double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
    202   double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
    203   void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
    204   bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration);
    205        
    206   bool CheckBounds(const Vector *x) const;
    207   void GetAlignvector(struct lsq_params * par) const;
    208 
    209   /// Initialising routines in fragmentation
    210   void CreateAdjacencyList2(ofstream *out, ifstream *output);
    211   void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
    212   void CreateListOfBondsPerAtom(ofstream *out);
    213 
    214   // Graph analysis
    215   MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack);
    216   void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
    217   bool PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack);
    218   bond * FindNextUnused(atom *vertex);
    219   void SetNextComponentNumber(atom *vertex, int nr);
    220   void InitComponentNumbers();
    221   void OutputComponentNumber(ofstream *out, atom *vertex);
    222   void ResetAllBondsToUnused();
    223   void ResetAllAtomNumbers();
    224   int CountCyclicBonds(ofstream *out);
    225   bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
    226   string GetColor(enum Shading color);
    227 
    228   molecule *CopyMolecule();
    229 
    230   /// Fragment molecule by two different approaches:
    231   int FragmentMolecule(ofstream *out, int Order, config *configuration);
    232   bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
    233   bool StoreAdjacencyToFile(ofstream *out, char *path);
    234   bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
    235   bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
    236   bool StoreOrderAtSiteFile(ofstream *out, char *path);
    237   bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
    238   bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
    239   bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
    240   bool CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex);
    241   bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
    242   void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
    243   /// -# BOSSANOVA
    244   void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
    245   int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
    246   bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
    247   molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
    248   void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
    249   int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
    250   int GuesstimateFragmentCount(ofstream *out, int order);
    251 
    252   // Recognize doubly appearing molecules in a list of them
    253   int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
    254   int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);
    255 
    256   // Output routines.
    257   bool Output(ofstream *out);
    258   bool OutputTrajectories(ofstream *out);
    259   void OutputListOfBonds(ofstream *out) const;
    260   bool OutputXYZ(ofstream *out) const;
    261   bool OutputTrajectoriesXYZ(ofstream *out);
    262   bool Checkout(ofstream *out) const;
    263   bool OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output);
    264 
    265   private:
    266   int last_atom;      //!< number given to last atom
    267   atom *InternalPointer;  //!< internal pointer for PointCloud
     204class molecule {
     205        public:
     206                double cell_size[6];//!< cell size
     207                periodentafel *elemente; //!< periodic table with each element
     208                atom *start;                            //!< start of atom list
     209                atom *end;                                      //!< end of atom list
     210                bond *first;                            //!< start of bond list
     211                bond *last;                              //!< end of bond list
     212                bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
     213                map<atom *, struct Trajectory> Trajectories; //!< contains old trajectory points
     214                int MDSteps;                            //!< The number of MD steps in Trajectories
     215                int *NumberOfBondsPerAtom;      //!< Number of Bonds each atom has
     216                int AtomCount;                                  //!< number of atoms, brought up-to-date by CountAtoms()
     217                int BondCount;                                  //!< number of atoms, brought up-to-date by CountBonds()
     218                int ElementCount;                        //!< how many unique elements are therein
     219                int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
     220                int NoNonHydrogen;      //!< number of non-hydrogen atoms in molecule
     221                int NoNonBonds;          //!< number of non-hydrogen bonds in molecule
     222                int NoCyclicBonds;      //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
     223                double BondDistance;    //!< typical bond distance used in CreateAdjacencyList() and furtheron
     224                bool ActiveFlag;    //!< in a MoleculeListClass used to discern active from inactive molecules
     225                Vector Center;      //!< Center of molecule in a global box
     226                char name[MAXSTRINGSIZE];         //!< arbitrary name
     227                int IndexNr;        //!< index of molecule in a MoleculeListClass
     228
     229        molecule(periodentafel *teil);
     230        ~molecule();
     231
     232        /// remove atoms from molecule.
     233        bool AddAtom(atom *pointer);
     234        bool RemoveAtom(atom *pointer);
     235        bool UnlinkAtom(atom *pointer);
     236        bool CleanupMolecule();
     237
     238        /// Add/remove atoms to/from molecule.
     239        atom * AddCopyAtom(atom *pointer);
     240        bool AddXYZFile(string filename);
     241        bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
     242        bond * AddBond(atom *first, atom *second, int degree);
     243        bool RemoveBond(bond *pointer);
     244        bool RemoveBonds(atom *BondPartner);
     245
     246        /// Find atoms.
     247        atom * FindAtom(int Nr) const;
     248        atom * AskAtom(string text);
     249
     250        /// Count and change present atoms' coordination.
     251        void CountAtoms(ofstream *out);
     252        void CountElements();
     253        void CalculateOrbitals(class config &configuration);
     254        bool CenterInBox(ofstream *out, Vector *BoxLengths);
     255        void CenterEdge(ofstream *out, Vector *max);
     256        void CenterOrigin(ofstream *out, Vector *max);
     257        void CenterGravity(ofstream *out, Vector *max);
     258        void Translate(const Vector *x);
     259        void Mirror(const Vector *x);
     260        void Align(Vector *n);
     261        void Scale(double **factor);
     262        void DetermineCenter(Vector &center);
     263        Vector * DetermineCenterOfGravity(ofstream *out);
     264        Vector * DetermineCenterOfAll(ofstream *out);
     265        void SetNameFromFilename(char *filename);
     266        void SetBoxDimension(Vector *dim);
     267        double * ReturnFullMatrixforSymmetric(double *cell_size);
     268        void ScanForPeriodicCorrection(ofstream *out);
     269        void PrincipalAxisSystem(ofstream *out, bool DoRotate);
     270        double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
     271        Vector* FindEmbeddingHole(ofstream *out, molecule *srcmol);
     272
     273        bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);
     274
     275        bool CheckBounds(const Vector *x) const;
     276        void GetAlignvector(struct lsq_params * par) const;
     277
     278        /// Initialising routines in fragmentation
     279        void CreateAdjacencyList2(ofstream *out, ifstream *output);
     280        void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
     281        void CreateListOfBondsPerAtom(ofstream *out);
     282
     283        // Graph analysis
     284        MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack);
     285        void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
     286        bool PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack);
     287        bond * FindNextUnused(atom *vertex);
     288        void SetNextComponentNumber(atom *vertex, int nr);
     289        void InitComponentNumbers();
     290        void OutputComponentNumber(ofstream *out, atom *vertex);
     291        void ResetAllBondsToUnused();
     292        void ResetAllAtomNumbers();
     293        int CountCyclicBonds(ofstream *out);
     294        bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
     295        string GetColor(enum Shading color);
     296
     297        molecule *CopyMolecule();
     298
     299        /// Fragment molecule by two different approaches:
     300        int FragmentMolecule(ofstream *out, int Order, config *configuration);
     301        bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
     302        bool StoreAdjacencyToFile(ofstream *out, char *path);
     303        bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
     304        bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
     305        bool StoreOrderAtSiteFile(ofstream *out, char *path);
     306        bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
     307        bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
     308        bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
     309        bool CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex);
     310        bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
     311        void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
     312        /// -# BOSSANOVA
     313        void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
     314        int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
     315        bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
     316        molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
     317        void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
     318        int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
     319        int GuesstimateFragmentCount(ofstream *out, int order);
     320
     321        // Recognize doubly appearing molecules in a list of them
     322        int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
     323        int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);
     324
     325        // Output routines.
     326        bool Output(ofstream *out);
     327        bool OutputTrajectories(ofstream *out);
     328        void OutputListOfBonds(ofstream *out) const;
     329        bool OutputXYZ(ofstream *out) const;
     330        bool OutputTrajectoriesXYZ(ofstream *out);
     331        bool Checkout(ofstream *out) const;
     332        bool OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output);
     333
     334        private:
     335        int last_atom;                  //!< number given to last atom
    268336};
    269337
     
    271339 */
    272340class MoleculeListClass {
    273   public:
    274     MoleculeList ListOfMolecules; //!< List of the contained molecules
    275     int MaxIndex;
    276 
    277   MoleculeListClass();
    278   ~MoleculeListClass();
    279 
    280   bool AddHydrogenCorrection(ofstream *out, char *path);
    281   bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
    282   void insert(molecule *mol);
    283   molecule * ReturnIndex(int index);
    284   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
    285   int NumberOfActiveMolecules();
    286   void Enumerate(ofstream *out);
    287   void Output(ofstream *out);
    288 
    289   // merging of molecules
     341        public:
     342          MoleculeList ListOfMolecules; //!< List of the contained molecules
     343          int MaxIndex;
     344
     345        MoleculeListClass();
     346        ~MoleculeListClass();
     347
     348        bool AddHydrogenCorrection(ofstream *out, char *path);
     349        bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
     350        bool insert(molecule *mol);
     351        molecule * ReturnIndex(int index);
     352        bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
     353        int NumberOfActiveMolecules();
     354        void Enumerate(ofstream *out);
     355        void Output(ofstream *out);
     356
     357        // merging of molecules
    290358  bool SimpleMerge(molecule *mol, molecule *srcmol);
    291359  bool SimpleAdd(molecule *mol, molecule *srcmol);
     
    295363  bool EmbedMerge(molecule *mol, molecule *srcmol);
    296364
    297   private:
     365        private:
    298366};
    299367
     
    303371 */
    304372class MoleculeLeafClass {
    305   public:
    306     molecule *Leaf;                   //!< molecule of this leaf
    307     //MoleculeLeafClass *UpLeaf;        //!< Leaf one level up
    308     //MoleculeLeafClass *DownLeaf;      //!< First leaf one level down
    309     MoleculeLeafClass *previous;  //!< Previous leaf on this level
    310     MoleculeLeafClass *next;      //!< Next leaf on this level
    311 
    312   //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
    313   MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
    314   ~MoleculeLeafClass();
    315 
    316   bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
    317   bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
    318   bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
    319   bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
    320   bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList);
    321   void TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
    322   int Count() const;
    323 };
    324 
    325 class ConfigFileBuffer {
    326   public:
    327     char **buffer;
    328     int *LineMapping;
    329     int CurrentLine;
    330     int NoLines;
    331 
    332     ConfigFileBuffer();
    333     ConfigFileBuffer(char *filename);
    334     ~ConfigFileBuffer();
    335 
    336     void InitMapping();
    337     void MapIonTypesInBuffer(int NoAtoms);
    338 };
    339 
     373        public:
     374                molecule *Leaf;                                                                  //!< molecule of this leaf
     375                //MoleculeLeafClass *UpLeaf;                            //!< Leaf one level up
     376                //MoleculeLeafClass *DownLeaf;                  //!< First leaf one level down
     377                MoleculeLeafClass *previous;    //!< Previous leaf on this level
     378                MoleculeLeafClass *next;                        //!< Next leaf on this level
     379
     380        //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
     381        MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
     382        ~MoleculeLeafClass();
     383
     384        bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
     385        bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
     386        bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
     387        bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
     388        bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList);
     389        void TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
     390        int Count() const;
     391};
     392
     393/** The config file.
     394 * The class contains all parameters that control a dft run also functions to load and save.
     395 */
     396class config {
     397        public:
     398                int PsiType;
     399                int MaxPsiDouble;
     400                int PsiMaxNoUp;
     401                int PsiMaxNoDown;
     402                int MaxMinStopStep;
     403                int InitMaxMinStopStep;
     404                int ProcPEGamma;
     405                int ProcPEPsi;
     406                char *configpath;
     407                char *configname;
     408                bool FastParsing;
     409                double Deltat;
     410                string basis;
     411
     412    char *databasepath;
     413
     414        private:
     415                char *mainname;
     416                char *defaultpath;
     417                char *pseudopotpath;
     418
     419                int DoOutVis;
     420                int DoOutMes;
     421                int DoOutNICS;
     422                int DoOutOrbitals;
     423                int DoOutCurrent;
     424                int DoFullCurrent;
     425                int DoPerturbation;
     426                int DoWannier;
     427                int CommonWannier;
     428                double SawtoothStart;
     429                int VectorPlane;
     430                double VectorCut;
     431                int UseAddGramSch;
     432                int Seed;
     433
     434                int MaxOuterStep;
     435                int OutVisStep;
     436                int OutSrcStep;
     437                double TargetTemp;
     438                int ScaleTempStep;
     439                int MaxPsiStep;
     440                double EpsWannier;
     441
     442                int MaxMinStep;
     443                double RelEpsTotalEnergy;
     444                double RelEpsKineticEnergy;
     445                int MaxMinGapStopStep;
     446                int MaxInitMinStep;
     447                double InitRelEpsTotalEnergy;
     448                double InitRelEpsKineticEnergy;
     449                int InitMaxMinGapStopStep;
     450
     451                //double BoxLength[NDIM*NDIM];
     452
     453                double ECut;
     454                int MaxLevel;
     455                int RiemannTensor;
     456                int LevRFactor;
     457                int RiemannLevel;
     458                int Lev0Factor;
     459                int RTActualUse;
     460                int AddPsis;
     461
     462                double RCut;
     463                int StructOpt;
     464                int IsAngstroem;
     465                int RelativeCoord;
     466                int MaxTypes;
     467
     468
     469        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);
     470
     471        public:
     472        config();
     473        ~config();
     474
     475        int TestSyntax(char *filename, periodentafel *periode, molecule *mol);
     476        void Load(char *filename, periodentafel *periode, molecule *mol);
     477        void LoadOld(char *filename, periodentafel *periode, molecule *mol);
     478        void RetrieveConfigPathAndName(string filename);
     479        bool Save(const char *filename, periodentafel *periode, molecule *mol) const;
     480        bool SaveMPQC(const char *filename, molecule *mol) const;
     481        void Edit();
     482        bool GetIsAngstroem() const;
     483        char *GetDefaultPath() const;
     484        void SetDefaultPath(const char *path);
     485};
    340486
    341487#endif /*MOLECULES_HPP_*/
Note: See TracChangeset for help on using the changeset viewer.