source: src/analysis_correlation.cpp@ 1da5f5

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 1da5f5 was 5108e1, checked in by Tillmann Crueger <crueger@…>, 15 years ago

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