source: molecuilder/src/molecules.hpp@ 451d7a

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

Gaussian basis for MPQC input files can now be specified with -B switch

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