source: src/analysis_correlation.cpp@ 88104f

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 88104f was 58bbd3, checked in by Frederik Heber <heber@…>, 15 years ago

Cases 'C' are now handled by CommandLineUI.

  • Property mode set to 100644
File size: 22.9 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 return outmap;
182};
183
184/** Calculates the distance (pair) correlation between a given element and a point.
185 * \param *molecules list of molecules structure
186 * \param &elements vector of elements to correlate with point
187 * \param *point vector to the correlation point
188 * \return Map of dobules with values as pairs of atom and the vector
189 */
190CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point )
191{
192 Info FunctionInfo(__func__);
193 CorrelationToPointMap *outmap = NULL;
194 double distance = 0.;
195 double *cell_size = World::getInstance().getDomain();
196
197 if (molecules->ListOfMolecules.empty()) {
198 DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
199 return outmap;
200 }
201 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
202 (*MolWalker)->doCountAtoms();
203 outmap = new CorrelationToPointMap;
204 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
205 if ((*MolWalker)->ActiveFlag) {
206 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
207 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
208 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
209 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
210 if ((*type == NULL) || ((*iter)->type == *type)) {
211 distance = (*iter)->node->PeriodicDistance(*point, cell_size);
212 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
213 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
214 }
215 }
216 }
217
218 return outmap;
219};
220
221/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
222 * \param *molecules list of molecules structure
223 * \param &elements vector of elements to correlate to point
224 * \param *point vector to the correlation point
225 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
226 * \return Map of dobules with values as pairs of atom and the vector
227 */
228CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
229{
230 Info FunctionInfo(__func__);
231 CorrelationToPointMap *outmap = NULL;
232 double distance = 0.;
233 int n[NDIM];
234 Vector periodicX;
235 Vector checkX;
236
237 if (molecules->ListOfMolecules.empty()) {
238 DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
239 return outmap;
240 }
241 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
242 (*MolWalker)->doCountAtoms();
243 outmap = new CorrelationToPointMap;
244 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
245 if ((*MolWalker)->ActiveFlag) {
246 double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
247 double * FullInverseMatrix = InverseMatrix(FullMatrix);
248 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
249 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
250 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
251 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
252 if ((*type == NULL) || ((*iter)->type == *type)) {
253 periodicX = *(*iter)->node;
254 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
255 // go through every range in xyz and get distance
256 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
257 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
258 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
259 checkX = Vector(n[0], n[1], n[2]) + periodicX;
260 checkX.MatrixMultiplication(FullMatrix);
261 distance = checkX.distance(*point);
262 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
263 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
264 }
265 }
266 }
267 delete[](FullMatrix);
268 delete[](FullInverseMatrix);
269 }
270
271 return outmap;
272};
273
274/** Calculates the distance (pair) correlation between a given element and a surface.
275 * \param *molecules list of molecules structure
276 * \param &elements vector of elements to correlate to surface
277 * \param *Surface pointer to Tesselation class surface
278 * \param *LC LinkedCell structure to quickly find neighbouring atoms
279 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
280 */
281CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
282{
283 Info FunctionInfo(__func__);
284 CorrelationToSurfaceMap *outmap = NULL;
285 double distance = 0;
286 class BoundaryTriangleSet *triangle = NULL;
287 Vector centroid;
288
289 if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
290 DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
291 return outmap;
292 }
293 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
294 (*MolWalker)->doCountAtoms();
295 outmap = new CorrelationToSurfaceMap;
296 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
297 if ((*MolWalker)->ActiveFlag) {
298 DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
299 if ((*MolWalker)->empty())
300 DoLog(1) && (1) && (Log() << Verbose(1) << "\t is empty." << endl);
301 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
302 DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
303 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
304 if ((*type == NULL) || ((*iter)->type == *type)) {
305 TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
306 distance = Intersections.GetSmallestDistance();
307 triangle = Intersections.GetClosestTriangle();
308 outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
309 }
310 }
311 } else {
312 DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
313 }
314
315 return outmap;
316};
317
318/** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
319 * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
320 * 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
321 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
322 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
323 * \param *molecules list of molecules structure
324 * \param &elements vector of elements to correlate to surface
325 * \param *Surface pointer to Tesselation class surface
326 * \param *LC LinkedCell structure to quickly find neighbouring atoms
327 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
328 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
329 */
330CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
331{
332 Info FunctionInfo(__func__);
333 CorrelationToSurfaceMap *outmap = NULL;
334 double distance = 0;
335 class BoundaryTriangleSet *triangle = NULL;
336 Vector centroid;
337 int n[NDIM];
338 Vector periodicX;
339 Vector checkX;
340
341 if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
342 DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
343 return outmap;
344 }
345 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
346 (*MolWalker)->doCountAtoms();
347 outmap = new CorrelationToSurfaceMap;
348 double ShortestDistance = 0.;
349 BoundaryTriangleSet *ShortestTriangle = NULL;
350 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
351 if ((*MolWalker)->ActiveFlag) {
352 double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
353 double * FullInverseMatrix = InverseMatrix(FullMatrix);
354 DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
355 for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
356 DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
357 for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
358 if ((*type == NULL) || ((*iter)->type == *type)) {
359 periodicX = *(*iter)->node;
360 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3
361 // go through every range in xyz and get distance
362 ShortestDistance = -1.;
363 for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
364 for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
365 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
366 checkX = Vector(n[0], n[1], n[2]) + periodicX;
367 checkX.MatrixMultiplication(FullMatrix);
368 TriangleIntersectionList Intersections(&checkX,Surface,LC);
369 distance = Intersections.GetSmallestDistance();
370 triangle = Intersections.GetClosestTriangle();
371 if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
372 ShortestDistance = distance;
373 ShortestTriangle = triangle;
374 }
375 }
376 // insert
377 outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
378 //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
379 }
380 }
381 delete[](FullMatrix);
382 delete[](FullInverseMatrix);
383 }
384
385 return outmap;
386};
387
388/** Returns the index of the bin for a given value.
389 * \param value value whose bin to look for
390 * \param BinWidth width of bin
391 * \param BinStart first bin
392 */
393int GetBin ( const double value, const double BinWidth, const double BinStart )
394{
395 Info FunctionInfo(__func__);
396 int bin =(int) (floor((value - BinStart)/BinWidth));
397 return (bin);
398};
399
400
401/** Prints correlation (double, int) pairs to file.
402 * \param *file file to write to
403 * \param *map map to write
404 */
405void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
406{
407 Info FunctionInfo(__func__);
408 *file << "BinStart\tCount" << endl;
409 for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
410 *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
411 }
412};
413
414/** Prints correlation (double, (atom*,atom*) ) pairs to file.
415 * \param *file file to write to
416 * \param *map map to write
417 */
418void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
419{
420 Info FunctionInfo(__func__);
421 *file << "BinStart\tAtom1\tAtom2" << endl;
422 for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
423 *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
424 }
425};
426
427/** Prints correlation (double, int) pairs to file.
428 * \param *file file to write to
429 * \param *map map to write
430 */
431void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
432{
433 Info FunctionInfo(__func__);
434 *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
435 for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
436 *file << runner->first;
437 for (int i=0;i<NDIM;i++)
438 *file << "\t" << setprecision(8) << (runner->second.first->node->at(i) - runner->second.second->at(i));
439 *file << endl;
440 }
441};
442
443/** Prints correlation (double, int) pairs to file.
444 * \param *file file to write to
445 * \param *map map to write
446 */
447void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
448{
449 Info FunctionInfo(__func__);
450 *file << "BinStart\tTriangle" << endl;
451 if (!map->empty())
452 for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
453 *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
454 }
455};
456
Note: See TracBrowser for help on using the repository browser.