source: src/analysis_correlation.cpp@ c78d44

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since c78d44 was c78d44, checked in by Frederik Heber <heber@…>, 15 years ago

Change to PairCorrelation...() functions, accept now vector<element *> instead of ptrs to element *.

  • modified ParseCommndLineOptions for all three cases 'C': Pair (E), Point (P) and Surface (S).
  • modified all Analysis(Pair)Correlation...() unit tests to have a vector of elements be delivered.
  • Property mode set to 100644
File size: 26.0 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 *molecules list of molecules structure
26 * \param &elements vector of elements to correlate
27 * \return Map of doubles with values the pair of the two atoms.
28 */
29PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements)
30{
31 Info FunctionInfo(__func__);
32 PairCorrelationMap *outmap = NULL;
33 double distance = 0.;
34 double *domain = World::getInstance().getDomain();
35
36 if (molecules->ListOfMolecules.empty()) {
37 DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
38 return outmap;
39 }
40 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
41 (*MolWalker)->doCountAtoms();
42
43 // create all possible pairs of elements
44 set <pair<element *, element *> > PairsOfElements;
45 if (elements.size() >= 2) {
46 for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
47 for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
48 if (type1 != type2) {
49 PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
50 DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
51 }
52 } else if (elements.size() == 1) { // one to all are valid
53 element *elemental = *elements.begin();
54 PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
55 PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
56 } else { // all elements valid
57 PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
58 }
59
60 outmap = new PairCorrelationMap;
61 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
62 if ((*MolWalker)->ActiveFlag) {
63 DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
64 eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
65 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
66 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
67 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
68 if ((*MolOtherWalker)->ActiveFlag) {
69 DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
70 for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
71 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
72 if ((*iter)->getId() < (*runner)->getId()){
73 for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
74 if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
75 distance = (*iter)->node->PeriodicDistance(*(*runner)->node, domain);
76 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
77 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
78 }
79 }
80 }
81 }
82 }
83 }
84 }
85 }
86 return outmap;
87};
88
89/** Calculates the pair correlation between given elements.
90 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
91 * \param *molecules list of molecules structure
92 * \param &elements vector of elements to correlate
93 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
94 * \return Map of doubles with values the pair of the two atoms.
95 */
96PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const int ranges[NDIM] )
97{
98 Info FunctionInfo(__func__);
99 PairCorrelationMap *outmap = NULL;
100 double distance = 0.;
101 int n[NDIM];
102 Vector checkX;
103 Vector periodicX;
104 int Othern[NDIM];
105 Vector checkOtherX;
106 Vector periodicOtherX;
107
108 if (molecules->ListOfMolecules.empty()) {
109 DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
110 return outmap;
111 }
112 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
113 (*MolWalker)->doCountAtoms();
114
115 // create all possible pairs of elements
116 set <pair<element *, element *> > PairsOfElements;
117 if (elements.size() >= 2) {
118 for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
119 for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
120 if (type1 != type2) {
121 PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
122 DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
123 }
124 } else if (elements.size() == 1) { // one to all are valid
125 element *elemental = *elements.begin();
126 PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
127 PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
128 } else { // all elements valid
129 PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
130 }
131
132 outmap = new PairCorrelationMap;
133 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
134 if ((*MolWalker)->ActiveFlag) {
135 double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
136 double * FullInverseMatrix = InverseMatrix(FullMatrix);
137 DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
138 eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
139 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
140 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
141 periodicX = *(*iter)->node;
142 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
143 // go through every range in xyz and get distance
144 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
145 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
146 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
147 checkX = Vector(n[0], n[1], n[2]) + periodicX;
148 checkX.MatrixMultiplication(FullMatrix);
149 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
150 if ((*MolOtherWalker)->ActiveFlag) {
151 DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
152 for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
153 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
154 if ((*iter)->getId() < (*runner)->getId()){
155 for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
156 if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
157 periodicOtherX = *(*runner)->node;
158 periodicOtherX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
159 // go through every range in xyz and get distance
160 for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
161 for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
162 for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
163 checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
164 checkOtherX.MatrixMultiplication(FullMatrix);
165 distance = checkX.distance(checkOtherX);
166 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
167 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
168 }
169 }
170 }
171 }
172 }
173 }
174 }
175 }
176 delete[](FullMatrix);
177 delete[](FullInverseMatrix);
178 }
179 }
180
181// outmap = new PairCorrelationMap;
182// for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
183// if ((*MolWalker)->ActiveFlag) {
184// double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
185// double * FullInverseMatrix = InverseMatrix(FullMatrix);
186// DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
187// for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
188// DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
189// if ((type1 == NULL) || ((*iter)->type == type1)) {
190// periodicX = *(*iter)->node;
191// periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
192// // go through every range in xyz and get distance
193// for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
194// for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
195// for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
196// checkX = Vector(n[0], n[1], n[2]) + periodicX;
197// checkX.MatrixMultiplication(FullMatrix);
198// for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
199// if ((*MolOtherWalker)->ActiveFlag) {
200// DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
201// for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
202// DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
203// if ((*iter)->nr < (*runner)->nr)
204// if ((type2 == NULL) || ((*runner)->type == type2)) {
205// periodicOtherX = *(*runner)->node;
206// periodicOtherX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
207// // go through every range in xyz and get distance
208// for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
209// for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
210// for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
211// checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
212// checkOtherX.MatrixMultiplication(FullMatrix);
213// distance = checkX.distance(checkOtherX);
214// //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
215// outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
216// }
217// }
218// }
219// }
220// }
221// }
222// }
223// delete[](FullMatrix);
224// delete[](FullInverseMatrix);
225// }
226
227 return outmap;
228};
229
230/** Calculates the distance (pair) correlation between a given element and a point.
231 * \param *molecules list of molecules structure
232 * \param &elements vector of elements to correlate with point
233 * \param *point vector to the correlation point
234 * \return Map of dobules with values as pairs of atom and the vector
235 */
236CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point )
237{
238 Info FunctionInfo(__func__);
239 CorrelationToPointMap *outmap = NULL;
240 double distance = 0.;
241
242 if (molecules->ListOfMolecules.empty()) {
243 DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
244 return outmap;
245 }
246 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
247 (*MolWalker)->doCountAtoms();
248 outmap = new CorrelationToPointMap;
249 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
250 if ((*MolWalker)->ActiveFlag) {
251 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
252 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
253 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
254 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
255 if ((*type == NULL) || ((*iter)->type == *type)) {
256 distance = (*iter)->node->PeriodicDistance(*point, World::getInstance().getDomain());
257 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
258 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
259 }
260 }
261 }
262
263 return outmap;
264};
265
266/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
267 * \param *molecules list of molecules structure
268 * \param &elements vector of elements to correlate to point
269 * \param *point vector to the correlation point
270 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
271 * \return Map of dobules with values as pairs of atom and the vector
272 */
273CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
274{
275 Info FunctionInfo(__func__);
276 CorrelationToPointMap *outmap = NULL;
277 double distance = 0.;
278 int n[NDIM];
279 Vector periodicX;
280 Vector checkX;
281
282 if (molecules->ListOfMolecules.empty()) {
283 DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
284 return outmap;
285 }
286 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
287 (*MolWalker)->doCountAtoms();
288 outmap = new CorrelationToPointMap;
289 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
290 if ((*MolWalker)->ActiveFlag) {
291 double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
292 double * FullInverseMatrix = InverseMatrix(FullMatrix);
293 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
294 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
295 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
296 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
297 if ((*type == NULL) || ((*iter)->type == *type)) {
298 periodicX = *(*iter)->node;
299 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
300 // go through every range in xyz and get distance
301 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
302 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
303 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
304 checkX = Vector(n[0], n[1], n[2]) + periodicX;
305 checkX.MatrixMultiplication(FullMatrix);
306 distance = checkX.distance(*point);
307 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
308 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
309 }
310 }
311 }
312 delete[](FullMatrix);
313 delete[](FullInverseMatrix);
314 }
315
316 return outmap;
317};
318
319/** Calculates the distance (pair) correlation between a given element and a surface.
320 * \param *molecules list of molecules structure
321 * \param &elements vector of elements to correlate to surface
322 * \param *Surface pointer to Tesselation class surface
323 * \param *LC LinkedCell structure to quickly find neighbouring atoms
324 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
325 */
326CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
327{
328 Info FunctionInfo(__func__);
329 CorrelationToSurfaceMap *outmap = NULL;
330 double distance = 0;
331 class BoundaryTriangleSet *triangle = NULL;
332 Vector centroid;
333
334 if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
335 DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
336 return outmap;
337 }
338 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
339 (*MolWalker)->doCountAtoms();
340 outmap = new CorrelationToSurfaceMap;
341 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
342 if ((*MolWalker)->ActiveFlag) {
343 DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
344 if ((*MolWalker)->empty())
345 DoLog(1) && (1) && (Log() << Verbose(1) << "\t is empty." << endl);
346 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
347 DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
348 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
349 if ((*type == NULL) || ((*iter)->type == *type)) {
350 TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
351 distance = Intersections.GetSmallestDistance();
352 triangle = Intersections.GetClosestTriangle();
353 outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
354 }
355 }
356 } else {
357 DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
358 }
359
360 return outmap;
361};
362
363/** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
364 * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
365 * 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
366 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
367 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
368 * \param *molecules list of molecules structure
369 * \param &elements vector of elements to correlate to surface
370 * \param *Surface pointer to Tesselation class surface
371 * \param *LC LinkedCell structure to quickly find neighbouring atoms
372 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
373 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
374 */
375CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
376{
377 Info FunctionInfo(__func__);
378 CorrelationToSurfaceMap *outmap = NULL;
379 double distance = 0;
380 class BoundaryTriangleSet *triangle = NULL;
381 Vector centroid;
382 int n[NDIM];
383 Vector periodicX;
384 Vector checkX;
385
386 if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
387 DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
388 return outmap;
389 }
390 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
391 (*MolWalker)->doCountAtoms();
392 outmap = new CorrelationToSurfaceMap;
393 double ShortestDistance = 0.;
394 BoundaryTriangleSet *ShortestTriangle = NULL;
395 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
396 if ((*MolWalker)->ActiveFlag) {
397 double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
398 double * FullInverseMatrix = InverseMatrix(FullMatrix);
399 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
400 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
401 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
402 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
403 if ((*type == NULL) || ((*iter)->type == *type)) {
404 periodicX = *(*iter)->node;
405 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
406 // go through every range in xyz and get distance
407 ShortestDistance = -1.;
408 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
409 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
410 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
411 checkX = Vector(n[0], n[1], n[2]) + periodicX;
412 checkX.MatrixMultiplication(FullMatrix);
413 TriangleIntersectionList Intersections(&checkX,Surface,LC);
414 distance = Intersections.GetSmallestDistance();
415 triangle = Intersections.GetClosestTriangle();
416 if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
417 ShortestDistance = distance;
418 ShortestTriangle = triangle;
419 }
420 }
421 // insert
422 outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
423 //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
424 }
425 }
426 delete[](FullMatrix);
427 delete[](FullInverseMatrix);
428 }
429
430 return outmap;
431};
432
433/** Returns the index of the bin for a given value.
434 * \param value value whose bin to look for
435 * \param BinWidth width of bin
436 * \param BinStart first bin
437 */
438int GetBin ( const double value, const double BinWidth, const double BinStart )
439{
440 Info FunctionInfo(__func__);
441 int bin =(int) (floor((value - BinStart)/BinWidth));
442 return (bin);
443};
444
445
446/** Prints correlation (double, int) pairs to file.
447 * \param *file file to write to
448 * \param *map map to write
449 */
450void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
451{
452 Info FunctionInfo(__func__);
453 *file << "BinStart\tCount" << endl;
454 for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
455 *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
456 }
457};
458
459/** Prints correlation (double, (atom*,atom*) ) pairs to file.
460 * \param *file file to write to
461 * \param *map map to write
462 */
463void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
464{
465 Info FunctionInfo(__func__);
466 *file << "BinStart\tAtom1\tAtom2" << endl;
467 for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
468 *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
469 }
470};
471
472/** Prints correlation (double, int) pairs to file.
473 * \param *file file to write to
474 * \param *map map to write
475 */
476void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
477{
478 Info FunctionInfo(__func__);
479 *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
480 for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
481 *file << runner->first;
482 for (int i=0;i<NDIM;i++)
483 *file << "\t" << setprecision(8) << (runner->second.first->node->at(i) - runner->second.second->at(i));
484 *file << endl;
485 }
486};
487
488/** Prints correlation (double, int) pairs to file.
489 * \param *file file to write to
490 * \param *map map to write
491 */
492void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
493{
494 Info FunctionInfo(__func__);
495 *file << "BinStart\tTriangle" << endl;
496 if (!map->empty())
497 for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
498 *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
499 }
500};
501
Note: See TracBrowser for help on using the repository browser.