source: molecuilder/src/molecule.hpp@ 7c2f6b

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

Complete refactoring of molecule_dynamics.cpp
  • Property mode set to 100755
File size: 19.3 KB
Line 
1/** \file molecule.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/*********************************************** includes ***********************************/
12
13// GSL headers
14#include <gsl/gsl_eigen.h>
15#include <gsl/gsl_heapsort.h>
16#include <gsl/gsl_linalg.h>
17#include <gsl/gsl_matrix.h>
18#include <gsl/gsl_multimin.h>
19#include <gsl/gsl_vector.h>
20#include <gsl/gsl_randist.h>
21
22//// STL headers
23#include <map>
24#include <set>
25#include <deque>
26#include <list>
27#include <vector>
28
29#include "graph.hpp"
30#include "stackclass.hpp"
31#include "tesselation.hpp"
32
33/****************************************** forward declarations *****************************/
34
35class atom;
36class bond;
37class element;
38class ForceMatrix;
39class LinkedCell;
40class molecule;
41class MoleculeLeafClass;
42class MoleculeListClass;
43class periodentafel;
44class Vector;
45
46/******************************** Some definitions for easier reading **********************************/
47
48#define MoleculeList list <molecule *>
49#define MoleculeListTest pair <MoleculeList::iterator, bool>
50
51#define DistancePair pair < double, atom* >
52#define DistanceMap multimap < double, atom* >
53#define DistanceTestPair pair < DistanceMap::iterator, bool>
54
55
56/************************************* Class definitions ****************************************/
57
58/** Structure to contain parameters needed for evaluation of constraint potential.
59 */
60struct EvaluatePotential
61{
62 int startstep; //!< start configuration (MDStep in atom::trajectory)
63 int endstep; //!< end configuration (MDStep in atom::trajectory)
64 atom **PermutationMap; //!< gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
65 DistanceMap **DistanceList; //!< distance list of each atom to each atom
66 DistanceMap::iterator *StepList; //!< iterator to ascend through NearestNeighbours \a **DistanceList
67 int *DoubleList; //!< count of which sources want to move to this target, basically the injective measure (>1 -> not injective)
68 DistanceMap::iterator *DistanceIterators; //!< marks which was the last picked target as injective candidate with smallest distance
69 bool IsAngstroem; //!< whether coordinates are in angstroem (true) or bohrradius (false)
70 double *PenaltyConstants; //!< penalty constant in front of each term
71};
72
73#define MaxThermostats 6 //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
74enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover }; //!< Thermostat names for output
75
76
77/** The complete molecule.
78 * Class incorporates number of types
79 */
80class molecule : public PointCloud {
81 public:
82 double cell_size[6];//!< cell size
83 periodentafel *elemente; //!< periodic table with each element
84 atom *start; //!< start of atom list
85 atom *end; //!< end of atom list
86 bond *first; //!< start of bond list
87 bond *last; //!< end of bond list
88 bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
89 int MDSteps; //!< The number of MD steps in Trajectories
90 int *NumberOfBondsPerAtom; //!< Number of Bonds each atom has
91 int AtomCount; //!< number of atoms, brought up-to-date by CountAtoms()
92 int BondCount; //!< number of atoms, brought up-to-date by CountBonds()
93 int ElementCount; //!< how many unique elements are therein
94 int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
95 int NoNonHydrogen; //!< number of non-hydrogen atoms in molecule
96 int NoNonBonds; //!< number of non-hydrogen bonds in molecule
97 int NoCyclicBonds; //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
98 double BondDistance; //!< typical bond distance used in CreateAdjacencyList() and furtheron
99 bool ActiveFlag; //!< in a MoleculeListClass used to discern active from inactive molecules
100 Vector Center; //!< Center of molecule in a global box
101 char name[MAXSTRINGSIZE]; //!< arbitrary name
102 int IndexNr; //!< index of molecule in a MoleculeListClass
103 class Tesselation *TesselStruct;
104
105 molecule(periodentafel *teil);
106 virtual ~molecule();
107
108 // re-definition of virtual functions from PointCloud
109 Vector *GetCenter(ofstream *out);
110 TesselPoint *GetPoint();
111 TesselPoint *GetTerminalPoint();
112 void GoToNext();
113 void GoToPrevious();
114 void GoToFirst();
115 void GoToLast();
116 bool IsEmpty();
117 bool IsEnd();
118
119 // templates for allowing global manipulation of all vectors
120 template <typename res> void ActOnAllVectors( res (Vector::*f)() );
121 template <typename res> void ActOnAllVectors( res (Vector::*f)() const);
122 template <typename res> void ActOnAllVectors( res (Vector::*f)() ) const;
123 template <typename res> void ActOnAllVectors( res (Vector::*f)() const) const;
124 template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t );
125 template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T) const, T t );
126 template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t ) const;
127 template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T) const, T t ) const;
128 template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u );
129 template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u );
130 template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const;
131 template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const;
132 template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v);
133 template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v);
134 template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const;
135 template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const;
136
137 // templates for allowing global manipulation of molecule with each atom as single argument
138 template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) );
139 template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) const);
140 template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) ) const;
141 template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) const) const;
142
143 // templates for allowing global copying of molecule with each atom as single argument
144 template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy);
145 template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const , molecule *copy);
146 template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const;
147 template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const, molecule *copy) const;
148
149 // templates for allowing global manipulation of all atoms
150 template <typename res> void ActOnAllAtoms( res (atom::*f)() );
151 template <typename res> void ActOnAllAtoms( res (atom::*f)() const );
152 template <typename res> void ActOnAllAtoms( res (atom::*f)() ) const;
153 template <typename res> void ActOnAllAtoms( res (atom::*f)() const) const;
154 template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T), T t );
155 template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T) const, T t );
156 template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T), T t ) const;
157 template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T) const, T t ) const;
158 template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U), T t, U u );
159 template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U) const, T t, U u );
160 template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U), T t, U u ) const;
161 template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U) const, T t, U u ) const;
162 template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v);
163 template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V) const, T t, U u, V v);
164 template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v) const;
165 template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V) const, T t, U u, V v) const;
166 template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W), T t, U u, V v, W w);
167 template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W) const, T t, U u, V v, W w);
168 template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W), T t, U u, V v, W w) const;
169 template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W) const, T t, U u, V v, W w) const;
170
171 // templates for allowing conditional global copying of molecule with each atom as single argument
172 template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () );
173 template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t );
174 template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u );
175 template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v );
176
177 // templates for allowing global manipulation of an array with one entry per atom
178 void SetIndexedArrayForEachAtomTo ( atom **array, int TesselPoint::* index);
179 template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, void (*Setor)(T *, T));
180 template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, void (*Setor)(T *, T), T t);
181 template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T));
182 template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T), T t);
183 template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, T (atom::*Setor)(Vector &), Vector atom::*value);
184 template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, T (atom::*Setor)(Vector &), Vector &vect );
185
186 /// remove atoms from molecule.
187 bool AddAtom(atom *pointer);
188 bool RemoveAtom(atom *pointer);
189 bool UnlinkAtom(atom *pointer);
190 bool CleanupMolecule();
191
192 /// Add/remove atoms to/from molecule.
193 atom * AddCopyAtom(atom *pointer);
194 bool AddXYZFile(string filename);
195 bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
196 bond * AddBond(atom *first, atom *second, int degree = 1);
197 bool RemoveBond(bond *pointer);
198 bool RemoveBonds(atom *BondPartner);
199
200 /// Find atoms.
201 atom * FindAtom(int Nr) const;
202 atom * AskAtom(string text);
203
204 /// Count and change present atoms' coordination.
205 void CountAtoms(ofstream *out);
206 void CountElements();
207 void CalculateOrbitals(class config &configuration);
208 bool CenterInBox(ofstream *out);
209 bool BoundInBox(ofstream *out);
210 void CenterEdge(ofstream *out, Vector *max);
211 void CenterOrigin(ofstream *out);
212 void CenterPeriodic(ofstream *out);
213 void CenterAtVector(ofstream *out, Vector *newcenter);
214 void Translate(const Vector *x);
215 void TranslatePeriodically(const Vector *trans);
216 void Mirror(const Vector *x);
217 void Align(Vector *n);
218 void Scale(double **factor);
219 void DeterminePeriodicCenter(Vector &center);
220 Vector * DetermineCenterOfGravity(ofstream *out);
221 Vector * DetermineCenterOfAll(ofstream *out);
222 void SetNameFromFilename(const char *filename);
223 void SetBoxDimension(Vector *dim);
224 void ScanForPeriodicCorrection(ofstream *out);
225 bool VerletForceIntegration(ofstream *out, char *file, config &configuration);
226 void Thermostats(config &configuration, double ActualTemp, int Thermostat);
227 void PrincipalAxisSystem(ofstream *out, bool DoRotate);
228 double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
229 Vector* FindEmbeddingHole(ofstream *out, molecule *srcmol);
230
231
232 double ConstrainedPotential(ofstream *out, struct EvaluatePotential &Params);
233 double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
234 void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
235 bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity);
236
237 bool CheckBounds(const Vector *x) const;
238 void GetAlignvector(struct lsq_params * par) const;
239
240 /// Initialising routines in fragmentation
241 void CreateAdjacencyList2(ofstream *out, ifstream *output);
242 void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
243 void CreateListOfBondsPerAtom(ofstream *out);
244
245 // Graph analysis
246 MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack);
247 void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
248 bool PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack);
249 bond * FindNextUnused(atom *vertex);
250 void SetNextComponentNumber(atom *vertex, int nr);
251 void InitComponentNumbers();
252 void OutputComponentNumber(ofstream *out, atom *vertex);
253 void ResetAllBondsToUnused();
254 void ResetAllAtomNumbers();
255 int CountCyclicBonds(ofstream *out);
256 bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
257 string GetColor(enum Shading color);
258
259 molecule *CopyMolecule();
260 molecule* CopyMoleculeFromSubRegion(Vector offset, double *parallelepiped);
261
262 /// Fragment molecule by two different approaches:
263 int FragmentMolecule(ofstream *out, int Order, config *configuration);
264 bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
265 bool StoreAdjacencyToFile(ofstream *out, char *path);
266 bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
267 bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
268 bool StoreOrderAtSiteFile(ofstream *out, char *path);
269 bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
270 bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
271 bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
272 bool CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex);
273 bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
274 void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
275 /// -# BOSSANOVA
276 void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
277 int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
278 bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
279 molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
280 void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
281 int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
282 int GuesstimateFragmentCount(ofstream *out, int order);
283
284 // Recognize doubly appearing molecules in a list of them
285 int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
286 int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);
287
288 // Output routines.
289 bool Output(ofstream *out);
290 bool OutputTrajectories(ofstream *out);
291 void OutputListOfBonds(ofstream *out) const;
292 bool OutputXYZ(ofstream *out) const;
293 bool OutputTrajectoriesXYZ(ofstream *out);
294 bool Checkout(ofstream *out) const;
295 bool OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output);
296
297 private:
298 int last_atom; //!< number given to last atom
299 atom *InternalPointer; //!< internal pointer for PointCloud
300};
301
302#include "molecule_template.hpp"
303
304/** A list of \a molecule classes.
305 */
306class MoleculeListClass {
307 public:
308 MoleculeList ListOfMolecules; //!< List of the contained molecules
309 int MaxIndex;
310
311 MoleculeListClass();
312 ~MoleculeListClass();
313
314 bool AddHydrogenCorrection(ofstream *out, char *path);
315 bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
316 void insert(molecule *mol);
317 molecule * ReturnIndex(int index);
318 bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
319 int NumberOfActiveMolecules();
320 void Enumerate(ofstream *out);
321 void Output(ofstream *out);
322
323 // merging of molecules
324 bool SimpleMerge(molecule *mol, molecule *srcmol);
325 bool SimpleAdd(molecule *mol, molecule *srcmol);
326 bool SimpleMultiMerge(molecule *mol, int *src, int N);
327 bool SimpleMultiAdd(molecule *mol, int *src, int N);
328 bool ScatterMerge(molecule *mol, int *src, int N);
329 bool EmbedMerge(molecule *mol, molecule *srcmol);
330
331 private:
332};
333
334
335/** A leaf for a tree of \a molecule class
336 * Wraps molecules in a tree structure
337 */
338class MoleculeLeafClass {
339 public:
340 molecule *Leaf; //!< molecule of this leaf
341 //MoleculeLeafClass *UpLeaf; //!< Leaf one level up
342 //MoleculeLeafClass *DownLeaf; //!< First leaf one level down
343 MoleculeLeafClass *previous; //!< Previous leaf on this level
344 MoleculeLeafClass *next; //!< Next leaf on this level
345
346 //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
347 MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
348 ~MoleculeLeafClass();
349
350 bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
351 bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
352 bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
353 bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
354 bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList);
355 void TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
356 int Count() const;
357};
358
359
360#endif /*MOLECULES_HPP_*/
361
Note: See TracBrowser for help on using the repository browser.