source: molecuilder/src/molecules.hpp@ a7b3b8

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

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