source: molecuilder/src/molecule.hpp@ 31ccb6

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

molecule::StoreAdjacencyToFile() and molecule::StoreBondsToFile() now take additional filename.

  • since ParseCommandLineOptions() has cases j and J which use the above functions, we have to generalize these functions to work also without a given path and with arbritrary filename.

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