source: molecuilder/src/molecules.hpp@ c830e8e

Last change on this file since c830e8e was c830e8e, checked in by Frederik Heber <heber@…>, 16 years ago

Basic implementation of Multiple molecules.

builder.cpp:

molecules.hpp:

moleculelist.cpp:

  • replaced listofmolecules array by STL list everywhere (only smaller changes necessary)
  • new merging function: SimpleMerge, SimpleAdd, SimpleMultiMerge, SimpleMultiAdd, (EmbedMerge, ScatterMerge ... both not finished). Add does not while merge does delete the src molecules.
  • new function: Enumerate(). Output of all molecules with number of atoms and chemical formula
  • new function: NumberOfActiveMolecules(). Counts the number of molecules in the list with ActiveFlag set.
  • new function: insert(). Inserts a molecule into the list with a unique index

molecules.cpp:

  • new function: SetNameFromFilename. Takes basename of a filename and sets name accordingly.
  • new function: UnlinkAtom. Only removes atom from list, does not delete it from memory.

atom.cpp:

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