source: molecuilder/src/analysis_correlation.cpp@ 1f2e46

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

Huge change: Log() << Verbose(.) --> DoLog(.) && (Log() << Verbose(.) << ...);

Most of the files are affected, but this is necessary as if DoLog() says verbosity is not enough, all the stream operators won"t get executed which saves substantial amount of computation time.

Signed-off-by: Frederik Heber <heber@…>
  • Property mode set to 100644
File size: 19.4 KB
Line 
1/*
2 * analysis.cpp
3 *
4 * Created on: Oct 13, 2009
5 * Author: heber
6 */
7
8#include <iostream>
9
10#include "analysis_correlation.hpp"
11#include "element.hpp"
12#include "info.hpp"
13#include "log.hpp"
14#include "molecule.hpp"
15#include "tesselation.hpp"
16#include "tesselationhelpers.hpp"
17#include "triangleintersectionlist.hpp"
18#include "vector.hpp"
19#include "verbose.hpp"
20#include "World.hpp"
21
22
23/** Calculates the pair correlation between given elements.
24 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
25 * \param *out output stream for debugging
26 * \param *molecules list of molecules structure
27 * \param *type1 first element or NULL (if any element)
28 * \param *type2 second element or NULL (if any element)
29 * \return Map of doubles with values the pair of the two atoms.
30 */
31PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
32{
33 Info FunctionInfo(__func__);
34 PairCorrelationMap *outmap = NULL;
35 double distance = 0.;
36
37 if (molecules->ListOfMolecules.empty()) {
38 DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
39 return outmap;
40 }
41 outmap = new PairCorrelationMap;
42 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
43 if ((*MolWalker)->ActiveFlag) {
44 DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
45 atom *Walker = (*MolWalker)->start;
46 while (Walker->next != (*MolWalker)->end) {
47 Walker = Walker->next;
48 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
49 if ((type1 == NULL) || (Walker->type == type1)) {
50 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
51 if ((*MolOtherWalker)->ActiveFlag) {
52 DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
53 atom *OtherWalker = (*MolOtherWalker)->start;
54 while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
55 OtherWalker = OtherWalker->next;
56 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
57 if (Walker->nr < OtherWalker->nr)
58 if ((type2 == NULL) || (OtherWalker->type == type2)) {
59 distance = Walker->node->PeriodicDistance(OtherWalker->node, World::get()->cell_size);
60 //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
61 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
62 }
63 }
64 }
65 }
66 }
67 }
68
69 return outmap;
70};
71
72/** Calculates the pair correlation between given elements.
73 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
74 * \param *out output stream for debugging
75 * \param *molecules list of molecules structure
76 * \param *type1 first element or NULL (if any element)
77 * \param *type2 second element or NULL (if any element)
78 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
79 * \return Map of doubles with values the pair of the two atoms.
80 */
81PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
82{
83 Info FunctionInfo(__func__);
84 PairCorrelationMap *outmap = NULL;
85 double distance = 0.;
86 int n[NDIM];
87 Vector checkX;
88 Vector periodicX;
89 int Othern[NDIM];
90 Vector checkOtherX;
91 Vector periodicOtherX;
92
93 if (molecules->ListOfMolecules.empty()) {
94 DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
95 return outmap;
96 }
97 outmap = new PairCorrelationMap;
98 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
99 if ((*MolWalker)->ActiveFlag) {
100 double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
101 double * FullInverseMatrix = InverseMatrix(FullMatrix);
102 DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
103 atom *Walker = (*MolWalker)->start;
104 while (Walker->next != (*MolWalker)->end) {
105 Walker = Walker->next;
106 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
107 if ((type1 == NULL) || (Walker->type == type1)) {
108 periodicX.CopyVector(Walker->node);
109 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
110 // go through every range in xyz and get distance
111 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
112 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
113 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
114 checkX.Init(n[0], n[1], n[2]);
115 checkX.AddVector(&periodicX);
116 checkX.MatrixMultiplication(FullMatrix);
117 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
118 if ((*MolOtherWalker)->ActiveFlag) {
119 DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
120 atom *OtherWalker = (*MolOtherWalker)->start;
121 while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
122 OtherWalker = OtherWalker->next;
123 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
124 if (Walker->nr < OtherWalker->nr)
125 if ((type2 == NULL) || (OtherWalker->type == type2)) {
126 periodicOtherX.CopyVector(OtherWalker->node);
127 periodicOtherX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
128 // go through every range in xyz and get distance
129 for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
130 for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
131 for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
132 checkOtherX.Init(Othern[0], Othern[1], Othern[2]);
133 checkOtherX.AddVector(&periodicOtherX);
134 checkOtherX.MatrixMultiplication(FullMatrix);
135 distance = checkX.Distance(&checkOtherX);
136 //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
137 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
138 }
139 }
140 }
141 }
142 }
143 }
144 }
145 Free(&FullMatrix);
146 Free(&FullInverseMatrix);
147 }
148
149 return outmap;
150};
151
152/** Calculates the distance (pair) correlation between a given element and a point.
153 * \param *out output stream for debugging
154 * \param *molecules list of molecules structure
155 * \param *type element or NULL (if any element)
156 * \param *point vector to the correlation point
157 * \return Map of dobules with values as pairs of atom and the vector
158 */
159CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
160{
161 Info FunctionInfo(__func__);
162 CorrelationToPointMap *outmap = NULL;
163 double distance = 0.;
164
165 if (molecules->ListOfMolecules.empty()) {
166 DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
167 return outmap;
168 }
169 outmap = new CorrelationToPointMap;
170 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
171 if ((*MolWalker)->ActiveFlag) {
172 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
173 atom *Walker = (*MolWalker)->start;
174 while (Walker->next != (*MolWalker)->end) {
175 Walker = Walker->next;
176 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
177 if ((type == NULL) || (Walker->type == type)) {
178 distance = Walker->node->PeriodicDistance(point, World::get()->cell_size);
179 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
180 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
181 }
182 }
183 }
184
185 return outmap;
186};
187
188/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
189 * \param *out output stream for debugging
190 * \param *molecules list of molecules structure
191 * \param *type element or NULL (if any element)
192 * \param *point vector to the correlation point
193 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
194 * \return Map of dobules with values as pairs of atom and the vector
195 */
196CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
197{
198 Info FunctionInfo(__func__);
199 CorrelationToPointMap *outmap = NULL;
200 double distance = 0.;
201 int n[NDIM];
202 Vector periodicX;
203 Vector checkX;
204
205 if (molecules->ListOfMolecules.empty()) {
206 DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
207 return outmap;
208 }
209 outmap = new CorrelationToPointMap;
210 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
211 if ((*MolWalker)->ActiveFlag) {
212 double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
213 double * FullInverseMatrix = InverseMatrix(FullMatrix);
214 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
215 atom *Walker = (*MolWalker)->start;
216 while (Walker->next != (*MolWalker)->end) {
217 Walker = Walker->next;
218 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
219 if ((type == NULL) || (Walker->type == type)) {
220 periodicX.CopyVector(Walker->node);
221 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
222 // go through every range in xyz and get distance
223 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
224 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
225 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
226 checkX.Init(n[0], n[1], n[2]);
227 checkX.AddVector(&periodicX);
228 checkX.MatrixMultiplication(FullMatrix);
229 distance = checkX.Distance(point);
230 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
231 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
232 }
233 }
234 }
235 Free(&FullMatrix);
236 Free(&FullInverseMatrix);
237 }
238
239 return outmap;
240};
241
242/** Calculates the distance (pair) correlation between a given element and a surface.
243 * \param *out output stream for debugging
244 * \param *molecules list of molecules structure
245 * \param *type element or NULL (if any element)
246 * \param *Surface pointer to Tesselation class surface
247 * \param *LC LinkedCell structure to quickly find neighbouring atoms
248 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
249 */
250CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
251{
252 Info FunctionInfo(__func__);
253 CorrelationToSurfaceMap *outmap = NULL;
254 double distance = 0;
255 class BoundaryTriangleSet *triangle = NULL;
256 Vector centroid;
257
258 if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
259 DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
260 return outmap;
261 }
262 outmap = new CorrelationToSurfaceMap;
263 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
264 if ((*MolWalker)->ActiveFlag) {
265 DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
266 atom *Walker = (*MolWalker)->start;
267 while (Walker->next != (*MolWalker)->end) {
268 Walker = Walker->next;
269 //Log() << Verbose(1) << "Current atom is " << *Walker << "." << endl;
270 if ((type == NULL) || (Walker->type == type)) {
271 TriangleIntersectionList Intersections(Walker->node,Surface,LC);
272 distance = Intersections.GetSmallestDistance();
273 triangle = Intersections.GetClosestTriangle();
274 outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
275 }
276 }
277 } else
278 DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
279
280
281 return outmap;
282};
283
284/** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
285 * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
286 * I.e. We multiply the atom::node with the inverse of the domain matrix, i.e. transform it to \f$[0,0^3\f$, then add per
287 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
288 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
289 * \param *out output stream for debugging
290 * \param *molecules list of molecules structure
291 * \param *type element or NULL (if any element)
292 * \param *Surface pointer to Tesselation class surface
293 * \param *LC LinkedCell structure to quickly find neighbouring atoms
294 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
295 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
296 */
297CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
298{
299 Info FunctionInfo(__func__);
300 CorrelationToSurfaceMap *outmap = NULL;
301 double distance = 0;
302 class BoundaryTriangleSet *triangle = NULL;
303 Vector centroid;
304 int n[NDIM];
305 Vector periodicX;
306 Vector checkX;
307
308 if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
309 DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
310 return outmap;
311 }
312 outmap = new CorrelationToSurfaceMap;
313 double ShortestDistance = 0.;
314 BoundaryTriangleSet *ShortestTriangle = NULL;
315 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
316 if ((*MolWalker)->ActiveFlag) {
317 double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
318 double * FullInverseMatrix = InverseMatrix(FullMatrix);
319 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
320 atom *Walker = (*MolWalker)->start;
321 while (Walker->next != (*MolWalker)->end) {
322 Walker = Walker->next;
323 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
324 if ((type == NULL) || (Walker->type == type)) {
325 periodicX.CopyVector(Walker->node);
326 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
327 // go through every range in xyz and get distance
328 ShortestDistance = -1.;
329 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
330 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
331 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
332 checkX.Init(n[0], n[1], n[2]);
333 checkX.AddVector(&periodicX);
334 checkX.MatrixMultiplication(FullMatrix);
335 TriangleIntersectionList Intersections(&checkX,Surface,LC);
336 distance = Intersections.GetSmallestDistance();
337 triangle = Intersections.GetClosestTriangle();
338 if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
339 ShortestDistance = distance;
340 ShortestTriangle = triangle;
341 }
342 }
343 // insert
344 outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
345 //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
346 }
347 }
348 Free(&FullMatrix);
349 Free(&FullInverseMatrix);
350 }
351
352 return outmap;
353};
354
355/** Returns the index of the bin for a given value.
356 * \param value value whose bin to look for
357 * \param BinWidth width of bin
358 * \param BinStart first bin
359 */
360int GetBin ( const double value, const double BinWidth, const double BinStart )
361{
362 Info FunctionInfo(__func__);
363 int bin =(int) (floor((value - BinStart)/BinWidth));
364 return (bin);
365};
366
367
368/** Prints correlation (double, int) pairs to file.
369 * \param *file file to write to
370 * \param *map map to write
371 */
372void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
373{
374 Info FunctionInfo(__func__);
375 *file << "BinStart\tCount" << endl;
376 for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
377 *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
378 }
379};
380
381/** Prints correlation (double, (atom*,atom*) ) pairs to file.
382 * \param *file file to write to
383 * \param *map map to write
384 */
385void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
386{
387 Info FunctionInfo(__func__);
388 *file << "BinStart\tAtom1\tAtom2" << endl;
389 for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
390 *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
391 }
392};
393
394/** Prints correlation (double, int) pairs to file.
395 * \param *file file to write to
396 * \param *map map to write
397 */
398void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
399{
400 Info FunctionInfo(__func__);
401 *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
402 for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
403 *file << runner->first;
404 for (int i=0;i<NDIM;i++)
405 *file << "\t" << setprecision(8) << (runner->second.first->node->x[i] - runner->second.second->x[i]);
406 *file << endl;
407 }
408};
409
410/** Prints correlation (double, int) pairs to file.
411 * \param *file file to write to
412 * \param *map map to write
413 */
414void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
415{
416 Info FunctionInfo(__func__);
417 *file << "BinStart\tTriangle" << endl;
418 if (!map->empty())
419 for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
420 *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
421 }
422};
423
Note: See TracBrowser for help on using the repository browser.