source: molecuilder/src/molecules.hpp@ 169b24

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

new functions MoleculeLeafClass::AssignKeySetsToFragment() and MoleculeLeafClass::FillListOfLocalAtoms() and molecule::OutputListOfBonds(), Count() now const
  • Property mode set to 100644
File size: 18.3 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_heapsort.h>
16
17// STL headers
18#include <map>
19#include <set>
20#include <deque>
21
22#include "helpers.hpp"
23#include "stackclass.hpp"
24
25class atom;
26class bond;
27class config;
28class element;
29class molecule;
30class MoleculeListClass;
31class periodentafel;
32class vector;
33class Verbose;
34
35/******************************** Some definitions for easier reading **********************************/
36
37#define KeyStack deque<int>
38#define KeySet set<int>
39#define Graph map<KeySet, pair<int, double>, KeyCompare >
40#define GraphPair pair<KeySet, pair<int, double> >
41#define KeySetTestPair pair<KeySet::iterator, bool>
42#define GraphTestPair pair<Graph::iterator, bool>
43
44struct KeyCompare
45{
46 bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
47};
48
49//bool operator < (KeySet SubgraphA, KeySet SubgraphB); //note: this declaration is important, otherwise normal < is used (producing wrong order)
50inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment); // Insert a KeySet into a Graph
51inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter); // Insert all KeySet's in a Graph into another Graph
52int CompareDoubles (const void * a, const void * b);
53
54
55/************************************* Class definitions ****************************************/
56
57/** Chemical element.
58 * Class incorporates data for a certain chemical element to be referenced from atom class.
59 */
60class element {
61 public:
62 double mass; //!< mass in g/mol
63 double CovalentRadius; //!< covalent radius
64 double VanDerWaalsRadius; //!< can-der-Waals radius
65 int Z; //!< atomic number
66 char name[64]; //!< atom name, i.e. "Hydrogren"
67 char symbol[3]; //!< short form of the atom, i.e. "H"
68 char period[8]; //!< period: n quantum number
69 char group[8]; //!< group: l quantum number
70 char block[8]; //!< block: l quantum number
71 element *previous; //!< previous item in list
72 element *next; //!< next element in list
73 int *sort; //!< sorc criteria
74 int No; //!< number of element set on periodentafel::Output()
75 double Valence; //!< number of valence electrons for this element
76 int NoValenceOrbitals; //!< number of valence orbitals, used for determining bond degree in molecule::CreateConnectmatrix()
77 double HBondDistance[NDIM]; //!< distance in Angstrom of this element to hydrogen (for single, double and triple bonds)
78 double HBondAngle[NDIM]; //!< typical angle for one, two, three bonded hydrogen (in degrees)
79
80 element();
81 ~element();
82
83 //> print element entries to screen
84 bool Output(ofstream *out) const;
85 bool Checkout(ofstream *out, const int No, const int NoOfAtoms) const;
86
87 private:
88};
89
90/** Periodentafel is a list of all elements sorted by their atomic number.
91 */
92class periodentafel {
93 public:
94 element *start; //!< start of element list
95 element *end; //!< end of element list
96 char header1[MAXSTRINGSIZE]; //!< store first header line
97 char header2[MAXSTRINGSIZE]; //!< store second header line
98
99 periodentafel();
100 ~periodentafel();
101
102 bool AddElement(element *pointer);
103 bool RemoveElement(element *pointer);
104 bool CleanupPeriodtable();
105 element * FindElement(int Z);
106 element * FindElement(char *shorthand) const;
107 element * AskElement();
108 bool Output(ofstream *output) const;
109 bool Checkout(ofstream *output, const int *checkliste) const;
110 bool LoadPeriodentafel(char *filename = NULL);
111 bool StorePeriodentafel(char *filename = NULL) const;
112
113 private:
114};
115
116// some algebraic matrix stuff
117#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
118#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2)) //!< hard-coded determinant of a 2x2 matrix
119
120/** Single vector.
121 * basically, just a x[3] but with helpful functions
122 */
123class vector {
124 public:
125 double x[NDIM];
126
127 vector();
128 ~vector();
129
130 double Distance(const vector *y) const;
131 double PeriodicDistance(const vector *y, const double *cell_size) const;
132 double ScalarProduct(const vector *y) const;
133 double Projection(const vector *y) const;
134 double Norm() const ;
135 double Angle(vector *y) const;
136
137 void AddVector(const vector *y);
138 void SubtractVector(const vector *y);
139 void CopyVector(const vector *y);
140 void RotateVector(const vector *y, const double alpha);
141 void Zero();
142 void Normalize();
143 void Translate(const vector *x);
144 void Mirror(const vector *x);
145 void Scale(double **factor);
146 void Scale(double *factor);
147 void Scale(double factor);
148 void MatrixMultiplication(double *M);
149 void InverseMatrixMultiplication(double *M);
150 void KeepPeriodic(ofstream *out, double *matrix);
151 void LinearCombinationOfVectors(const vector *x1, const vector *x2, const vector *x3, double *factors);
152
153 bool GetOneNormalVector(const vector *x1);
154 bool MakeNormalVector(const vector *y1);
155 bool MakeNormalVector(const vector *y1, const vector *y2);
156 bool MakeNormalVector(const vector *x1, const vector *x2, const vector *x3);
157 bool SolveSystem(vector *x1, vector *x2, vector *y, double alpha, double beta, double c);
158 bool LSQdistance(vector **vectors, int dim);
159
160 void AskPosition(double *cell_size, bool check);
161 bool Output(ofstream *out) const;
162};
163
164ofstream& operator<<(ofstream& ost, vector& m);
165
166/** Parameter structure for least square minimsation.
167 */
168struct LSQ_params {
169 vector **vectors;
170 int num;
171};
172
173double LSQ(const gsl_vector * x, void * params);
174
175/** Parameter structure for least square minimsation.
176 */
177struct lsq_params {
178 gsl_vector *x;
179 const molecule *mol;
180 element *type;
181};
182
183
184
185/** Single atom.
186 * Class incoporates position, type
187 */
188class atom {
189 public:
190 vector x; //!< coordinate array of atom, giving position within cell
191 vector v; //!< velocity array of atom
192 element *type; //!< pointing to element
193 atom *previous; //!< previous atom in molecule list
194 atom *next; //!< next atom in molecule list
195 atom *father; //!< In many-body bond order fragmentations points to originating atom
196 atom *Ancestor; //!< "Father" in Depth-First-Search
197 char *Name; //!< unique name used during many-body bond-order fragmentation
198 int FixedIon; //!< config variable that states whether forces act on the ion or not
199 int *sort; //!< sort criteria
200 int nr; //!< continuous, unique number
201 int GraphNr; //!< unique number, given in DepthFirstSearchAnalysis()
202 int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
203 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.
204 bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, given in DepthFirstSearchAnalysis()
205 unsigned char AdaptiveOrder; //!< current present bond order at site (0 means "not set")
206
207 atom();
208 ~atom();
209
210 bool Output(int ElementNo, int AtomNo, ofstream *out) const;
211 bool OutputXYZLine(ofstream *out) const;
212 atom *GetTrueFather();
213 bool Compare(atom &ptr);
214
215 private:
216};
217
218ostream & operator << (ostream &ost, atom &a);
219
220/** Bonds between atoms.
221 * Class incorporates bonds between atoms in a molecule,
222 * used to derive tge fragments in many-body bond order
223 * calculations.
224 */
225class bond {
226 public:
227 atom *leftatom; //!< first bond partner
228 atom *rightatom; //!< second bond partner
229 bond *previous; //!< previous atom in molecule list
230 bond *next; //!< next atom in molecule list
231 int HydrogenBond; //!< Number of hydrogen atoms in the bond
232 int BondDegree; //!< single, double, triple, ... bond
233 int nr; //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
234 bool Cyclic; //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
235 enum EdgeType Type;//!< whether this is a tree or back edge
236
237 atom * GetOtherAtom(atom *Atom) const;
238 bond * GetFirstBond();
239 bond * GetLastBond();
240
241 bool MarkUsed(enum Shading color);
242 enum Shading IsUsed();
243 void ResetUsed();
244 bool Contains(const atom *ptr);
245 bool Contains(const int nr);
246
247 bond();
248 bond(atom *left, atom *right);
249 bond(atom *left, atom *right, int degree);
250 bond(atom *left, atom *right, int degree, int number);
251 ~bond();
252
253 private:
254 enum Shading Used; //!< marker in depth-first search, DepthFirstSearchAnalysis()
255};
256
257ostream & operator << (ostream &ost, bond &b);
258
259class MoleculeLeafClass;
260
261/** The complete molecule.
262 * Class incorporates number of types
263 */
264class molecule {
265 public:
266 double cell_size[6];//!< cell size
267 periodentafel *elemente; //!< periodic table with each element
268 atom *start; //!< start of atom list
269 atom *end; //!< end of atom list
270 bond *first; //!< start of bond list
271 bond *last; //!< end of bond list
272 bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
273 int *NumberOfBondsPerAtom; //!< Number of Bonds each atom has
274 int AtomCount; //!< number of atoms, brought up-to-date by CountAtoms()
275 int BondCount; //!< number of atoms, brought up-to-date by CountBonds()
276 int ElementCount; //!< how many unique elements are therein
277 int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
278 int NoNonHydrogen; //!< number of non-hydrogen atoms in molecule
279 int NoNonBonds; //!< number of non-hydrogen bonds in molecule
280 int NoCyclicBonds; //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
281 double BondDistance; //!< typical bond distance used in CreateAdjacencyList() and furtheron
282
283 molecule(periodentafel *teil);
284 ~molecule();
285
286 /// remove atoms from molecule.
287 bool AddAtom(atom *pointer);
288 bool RemoveAtom(atom *pointer);
289 bool CleanupMolecule();
290
291 /// Add/remove atoms to/from molecule.
292 atom * AddCopyAtom(atom *pointer);
293 bool AddXYZFile(string filename);
294 bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
295 bond * AddBond(atom *first, atom *second, int degree);
296 bool RemoveBond(bond *pointer);
297 bool RemoveBonds(atom *BondPartner);
298
299 /// Find atoms.
300 atom * FindAtom(int Nr) const;
301 atom * AskAtom(char *text);
302
303 /// Count and change present atoms' coordination.
304 void CountAtoms(ofstream *out);
305 void CountElements();
306 void CalculateOrbitals(class config &configuration);
307 void CenterEdge(ofstream *out, vector *max);
308 void CenterOrigin(ofstream *out, vector *max);
309 void CenterGravity(ofstream *out, vector *max);
310 void Translate(const vector *x);
311 void Mirror(const vector *x);
312 void Align(vector *n);
313 void Scale(double **factor);
314 void DetermineCenterOfGravity(vector &CenterOfGravity);
315 void SetBoxDimension(vector *dim);
316 double * ReturnFullMatrixforSymmetric(double *cell_size);
317 void ScanForPeriodicCorrection(ofstream *out);
318
319 bool CheckBounds(const vector *x) const;
320 void GetAlignVector(struct lsq_params * par) const;
321
322 /// Initialising routines in fragmentation
323 void CreateAdjacencyList(ofstream *out, double bonddistance);
324 void CreateListOfBondsPerAtom(ofstream *out);
325
326 // Graph analysis
327 MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int &MinimumRingSize);
328 void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int &MinimumRingSize);
329 bond * FindNextUnused(atom *vertex);
330 void SetNextComponentNumber(atom *vertex, int nr);
331 void InitComponentNumbers();
332 void OutputComponentNumber(ofstream *out, atom *vertex);
333 void ResetAllBondsToUnused();
334 void ResetAllAtomNumbers();
335 int CountCyclicBonds(ofstream *out);
336 char * GetColor(enum Shading color);
337
338 molecule *CopyMolecule();
339
340 /// Fragment molecule by two different approaches:
341 void FragmentMolecule(ofstream *out, int Order, config *configuration);
342 void CheckOrderAtSite(ofstream *out, bool &FragmentationToDo, int Order);
343 bool StoreAdjacencyToFile(ofstream *out, char *path);
344 bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
345 bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
346 bool StoreOrderAtSiteFile(ofstream *out, char *path);
347 bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList, bool IsAngstroem);
348 bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
349 bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
350 bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
351 void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
352 /// -# BOSSANOVA
353 void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack);
354 int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
355 bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
356 molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
357 void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
358 int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
359 int GuesstimateFragmentCount(ofstream *out, int order);
360
361 // Recognize doubly appearing molecules in a list of them
362 int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
363 int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);
364
365 // Output routines.
366 bool Output(ofstream *out);
367 void OutputListOfBonds(ofstream *out) const;
368 bool OutputXYZ(ofstream *out) const;
369 bool Checkout(ofstream *out) const;
370
371 private:
372 int last_atom; //!< number given to last atom
373};
374
375/** A list of \a molecule classes.
376 */
377class MoleculeListClass {
378 public:
379 molecule **ListOfMolecules; //!< pointer list of fragment molecules to check for equality
380 int NumberOfMolecules; //!< Number of entries in \a **FragmentList and of to be returned one.
381 int NumberOfTopAtoms; //!< Number of atoms in the molecule from which all fragments originate
382
383 MoleculeListClass();
384 MoleculeListClass(int Num, int NumAtoms);
385 ~MoleculeListClass();
386
387 /// Output configs.
388 bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
389 bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
390 void Output(ofstream *out);
391
392 private:
393};
394
395
396/** A leaf for a tree of \a molecule class
397 * Wraps molecules in a tree structure
398 */
399class MoleculeLeafClass {
400 public:
401 molecule *Leaf; //!< molecule of this leaf
402 //MoleculeLeafClass *UpLeaf; //!< Leaf one level up
403 //MoleculeLeafClass *DownLeaf; //!< First leaf one level down
404 MoleculeLeafClass *previous; //!< Previous leaf on this level
405 MoleculeLeafClass *next; //!< Next leaf on this level
406
407 //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
408 MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
409 ~MoleculeLeafClass();
410
411 bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
412 bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
413 bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, int Order, int &FragmentCounter);
414 bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
415 bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, int &FragmentCounter, int GlobalAtomCount, bool &FreeList);
416 int Count() const;
417};
418
419/** The config file.
420 * The class contains all parameters that control a dft run also functions to load and save.
421 */
422class config {
423 public:
424 int PsiType;
425 int MaxPsiDouble;
426 int PsiMaxNoUp;
427 int PsiMaxNoDown;
428 int MaxMinStopStep;
429 int InitMaxMinStopStep;
430 int ProcPEGamma;
431 int ProcPEPsi;
432 char *configpath;
433 char *configname;
434
435 private:
436 char *mainname;
437 char *defaultpath;
438 char *pseudopotpath;
439
440 int DoOutVis;
441 int DoOutMes;
442 int DoOutNICS;
443 int DoOutOrbitals;
444 int DoOutCurrent;
445 int DoFullCurrent;
446 int DoPerturbation;
447 int CommonWannier;
448 double SawtoothStart;
449 int VectorPlane;
450 double VectorCut;
451 int UseAddGramSch;
452 int Seed;
453
454 int MaxOuterStep;
455 double Deltat;
456 int OutVisStep;
457 int OutSrcStep;
458 double TargetTemp;
459 int ScaleTempStep;
460 int MaxPsiStep;
461 double EpsWannier;
462
463 int MaxMinStep;
464 double RelEpsTotalEnergy;
465 double RelEpsKineticEnergy;
466 int MaxMinGapStopStep;
467 int MaxInitMinStep;
468 double InitRelEpsTotalEnergy;
469 double InitRelEpsKineticEnergy;
470 int InitMaxMinGapStopStep;
471
472 //double BoxLength[NDIM*NDIM];
473
474 double ECut;
475 int MaxLevel;
476 int RiemannTensor;
477 int LevRFactor;
478 int RiemannLevel;
479 int Lev0Factor;
480 int RTActualUse;
481 int AddPsis;
482
483 double RCut;
484 int StructOpt;
485 int IsAngstroem;
486 int RelativeCoord;
487 int MaxTypes;
488
489
490 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);
491
492 public:
493 config();
494 ~config();
495
496 int TestSyntax(char *filename, periodentafel *periode, molecule *mol);
497 void Load(char *filename, periodentafel *periode, molecule *mol);
498 void LoadOld(char *filename, periodentafel *periode, molecule *mol);
499 void RetrieveConfigPathAndName(char *filename);
500 bool Save(ofstream *file, periodentafel *periode, molecule *mol) const;
501 void Edit(molecule *mol);
502 bool GetIsAngstroem() const;
503 char *GetDefaultPath() const;
504 void config::SetDefaultPath(const char *path);
505};
506
507#endif /*MOLECULES_HPP_*/
508
Note: See TracBrowser for help on using the repository browser.