source: src/linkedcell.cpp@ ba9f5b

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 ba9f5b was bdc91e, checked in by Frederik Heber <heber@…>, 14 years ago

MEMFIXES: Tesselation routines were leaking memory.

Signed-off-by: Frederik Heber <heber@…>

  • Property mode set to 100644
File size: 13.9 KB
Line 
1/** \file linkedcell.cpp
2 *
3 * Function implementations for the class LinkedCell.
4 *
5 */
6
7#include "Helpers/MemDebug.hpp"
8
9#include "atom.hpp"
10#include "helpers.hpp"
11#include "linkedcell.hpp"
12#include "log.hpp"
13#include "molecule.hpp"
14#include "tesselation.hpp"
15#include "vector.hpp"
16
17// ========================================================= class LinkedCell ===========================================
18
19
20/** Constructor for class LinkedCell.
21 */
22LinkedCell::LinkedCell()
23{
24 LC = NULL;
25 for(int i=0;i<NDIM;i++)
26 N[i] = 0;
27 index = -1;
28 RADIUS = 0.;
29 max.Zero();
30 min.Zero();
31};
32
33/** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
34 * \param *set LCNodeSet class with all LCNode's
35 * \param RADIUS edge length of cells
36 */
37LinkedCell::LinkedCell(const PointCloud * const set, const double radius)
38{
39 TesselPoint *Walker = NULL;
40
41 RADIUS = radius;
42 LC = NULL;
43 for(int i=0;i<NDIM;i++)
44 N[i] = 0;
45 index = -1;
46 max.Zero();
47 min.Zero();
48 DoLog(1) && (Log() << Verbose(1) << "Begin of LinkedCell" << endl);
49 if ((set == NULL) || (set->IsEmpty())) {
50 DoeLog(1) && (eLog()<< Verbose(1) << "set is NULL or contains no linked cell nodes!" << endl);
51 return;
52 }
53 // 1. find max and min per axis of atoms
54 set->GoToFirst();
55 Walker = set->GetPoint();
56 for (int i=0;i<NDIM;i++) {
57 max[i] = Walker->node->at(i);
58 min[i] = Walker->node->at(i);
59 }
60 set->GoToFirst();
61 while (!set->IsEnd()) {
62 Walker = set->GetPoint();
63 for (int i=0;i<NDIM;i++) {
64 if (max[i] < Walker->node->at(i))
65 max[i] = Walker->node->at(i);
66 if (min[i] > Walker->node->at(i))
67 min[i] = Walker->node->at(i);
68 }
69 set->GoToNext();
70 }
71 DoLog(2) && (Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl);
72
73 // 2. find then number of cells per axis
74 for (int i=0;i<NDIM;i++) {
75 N[i] = static_cast<int>(floor((max[i] - min[i])/RADIUS)+1);
76 }
77 DoLog(2) && (Log() << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl);
78
79 // 3. allocate the lists
80 DoLog(2) && (Log() << Verbose(2) << "Allocating cells ... ");
81 if (LC != NULL) {
82 DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl);
83 return;
84 }
85 LC = new LinkedNodes[N[0]*N[1]*N[2]];
86 for (index=0;index<N[0]*N[1]*N[2];index++) {
87 LC [index].clear();
88 }
89 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
90
91 // 4. put each atom into its respective cell
92 DoLog(2) && (Log() << Verbose(2) << "Filling cells ... ");
93 set->GoToFirst();
94 while (!set->IsEnd()) {
95 Walker = set->GetPoint();
96 for (int i=0;i<NDIM;i++) {
97 n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
98 }
99 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
100 LC[index].push_back(Walker);
101 //Log() << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
102 set->GoToNext();
103 }
104 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
105 DoLog(1) && (Log() << Verbose(1) << "End of LinkedCell" << endl);
106};
107
108
109/** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
110 * \param *set LCNodeSet class with all LCNode's
111 * \param RADIUS edge length of cells
112 */
113LinkedCell::LinkedCell(LinkedNodes *set, const double radius)
114{
115 class TesselPoint *Walker = NULL;
116 RADIUS = radius;
117 LC = NULL;
118 for(int i=0;i<NDIM;i++)
119 N[i] = 0;
120 index = -1;
121 max.Zero();
122 min.Zero();
123 DoLog(1) && (Log() << Verbose(1) << "Begin of LinkedCell" << endl);
124 if (set->empty()) {
125 DoeLog(1) && (eLog()<< Verbose(1) << "set contains no linked cell nodes!" << endl);
126 return;
127 }
128 // 1. find max and min per axis of atoms
129 LinkedNodes::iterator Runner = set->begin();
130 for (int i=0;i<NDIM;i++) {
131 max[i] = (*Runner)->node->at(i);
132 min[i] = (*Runner)->node->at(i);
133 }
134 for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
135 Walker = *Runner;
136 for (int i=0;i<NDIM;i++) {
137 if (max[i] < Walker->node->at(i))
138 max[i] = Walker->node->at(i);
139 if (min[i] > Walker->node->at(i))
140 min[i] = Walker->node->at(i);
141 }
142 }
143 DoLog(2) && (Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl);
144
145 // 2. find then number of cells per axis
146 for (int i=0;i<NDIM;i++) {
147 N[i] = static_cast<int>(floor((max[i] - min[i])/RADIUS)+1);
148 }
149 DoLog(2) && (Log() << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl);
150
151 // 3. allocate the lists
152 DoLog(2) && (Log() << Verbose(2) << "Allocating cells ... ");
153 if (LC != NULL) {
154 DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl);
155 return;
156 }
157 LC = new LinkedNodes[N[0]*N[1]*N[2]];
158 for (index=0;index<N[0]*N[1]*N[2];index++) {
159 LC [index].clear();
160 }
161 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
162
163 // 4. put each atom into its respective cell
164 DoLog(2) && (Log() << Verbose(2) << "Filling cells ... ");
165 for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
166 Walker = *Runner;
167 for (int i=0;i<NDIM;i++) {
168 n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
169 }
170 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
171 LC[index].push_back(Walker);
172 //Log() << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
173 }
174 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
175 DoLog(1) && (Log() << Verbose(1) << "End of LinkedCell" << endl);
176};
177
178/** Destructor for class LinkedCell.
179 */
180LinkedCell::~LinkedCell()
181{
182 if (LC != NULL)
183 for (index=0;index<N[0]*N[1]*N[2];index++)
184 LC[index].clear();
185 delete[](LC);
186 for(int i=0;i<NDIM;i++)
187 N[i] = 0;
188 index = -1;
189};
190
191/** Checks whether LinkedCell::n[] is each within [0,N[]].
192 * \return if all in intervals - true, else -false
193 */
194bool LinkedCell::CheckBounds() const
195{
196 bool status = true;
197 for(int i=0;i<NDIM;i++)
198 status = status && ((n[i] >=0) && (n[i] < N[i]));
199// if (!status)
200// DoeLog(1) && (eLog()<< Verbose(1) << "indices are out of bounds!" << endl);
201 return status;
202};
203
204/** Checks whether LinkedCell::n[] plus relative offset is each within [0,N[]].
205 * Note that for this check we don't admonish if out of bounds.
206 * \param relative[NDIM] relative offset to current cell
207 * \return if all in intervals - true, else -false
208 */
209bool LinkedCell::CheckBounds(const int relative[NDIM]) const
210{
211 bool status = true;
212 for(int i=0;i<NDIM;i++)
213 status = status && ((n[i]+relative[i] >=0) && (n[i]+relative[i] < N[i]));
214 return status;
215};
216
217
218/** Returns a pointer to the current cell.
219 * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
220 */
221const LinkedCell::LinkedNodes* LinkedCell::GetCurrentCell() const
222{
223 if (CheckBounds()) {
224 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
225 return (&(LC[index]));
226 } else {
227 return NULL;
228 }
229};
230
231/** Returns a pointer to the current cell.
232 * \param relative[NDIM] offset for each axis with respect to the current cell LinkedCell::n[NDIM]
233 * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[]+relative[] are out of bounds.
234 */
235const LinkedCell::LinkedNodes* LinkedCell::GetRelativeToCurrentCell(const int relative[NDIM]) const
236{
237 if (CheckBounds(relative)) {
238 index = (n[0]+relative[0]) * N[1] * N[2] + (n[1]+relative[1]) * N[2] + (n[2]+relative[2]);
239 return (&(LC[index]));
240 } else {
241 return NULL;
242 }
243};
244
245/** Set the index to the cell containing a given Vector *x.
246 * \param *x Vector with coordinates
247 * \return Vector is inside bounding box - true, else - false
248 */
249bool LinkedCell::SetIndexToVector(const Vector * const x) const
250{
251 for (int i=0;i<NDIM;i++)
252 n[i] = (int)floor((x->at(i) - min[i])/RADIUS);
253
254 return CheckBounds();
255};
256
257/** Calculates the index for a given LCNode *Walker.
258 * \param *Walker LCNode to set index tos
259 * \return if the atom is also found in this cell - true, else - false
260 */
261bool LinkedCell::SetIndexToNode(const TesselPoint * const Walker) const
262{
263 bool status = false;
264 for (int i=0;i<NDIM;i++) {
265 n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
266 }
267 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
268 if (CheckBounds()) {
269 for (LinkedNodes::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
270 status = status || ((*Runner) == Walker);
271 return status;
272 } else {
273 DoeLog(1) && (eLog()<< Verbose(1) << "Node at " << *Walker << " is out of bounds." << endl);
274 return false;
275 }
276};
277
278/** Calculates the interval bounds of the linked cell grid.
279 * \param lower lower bounds
280 * \param upper upper bounds
281 * \param step how deep to check the neighbouring cells (i.e. number of layers to check)
282 */
283void LinkedCell::GetNeighbourBounds(int lower[NDIM], int upper[NDIM], int step) const
284{
285 for (int i=0;i<NDIM;i++) {
286 lower[i] = n[i]-step;
287 if (lower[i] < 0)
288 lower[i] = 0;
289 if (lower[i] >= N[i])
290 lower[i] = N[i]-1;
291 upper[i] = n[i]+step;
292 if (upper[i] >= N[i])
293 upper[i] = N[i]-1;
294 if (upper[i] < 0)
295 upper[i] = 0;
296 //Log() << Verbose(0) << "axis " << i << " has bounds [" << lower[i] << "," << upper[i] << "]" << endl;
297 }
298};
299
300/** Returns a list with all neighbours from the current LinkedCell::index.
301 * \param distance (if no distance, then adjacent cells are taken)
302 * \return list of tesselpoints
303 */
304LinkedCell::LinkedNodes* LinkedCell::GetallNeighbours(const double distance) const
305{
306 int Nlower[NDIM], Nupper[NDIM];
307 TesselPoint *Walker = NULL;
308 LinkedNodes *TesselList = new LinkedNodes;
309
310 // then go through the current and all neighbouring cells and check the contained points for possible candidates
311 const int step = (distance == 0) ? 1 : (int)floor(distance/RADIUS + 1.);
312 GetNeighbourBounds(Nlower, Nupper, step);
313
314 //Log() << Verbose(0) << endl;
315 for (n[0] = Nlower[0]; n[0] <= Nupper[0]; n[0]++)
316 for (n[1] = Nlower[1]; n[1] <= Nupper[1]; n[1]++)
317 for (n[2] = Nlower[2]; n[2] <= Nupper[2]; n[2]++) {
318 const LinkedNodes *List = GetCurrentCell();
319 //Log() << Verbose(1) << "Current cell is " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
320 if (List != NULL) {
321 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
322 Walker = *Runner;
323 TesselList->push_back(Walker);
324 }
325 }
326 }
327 return TesselList;
328};
329
330/** Set the index to the cell containing a given Vector *x, which is not inside the LinkedCell's domain
331 * Note that as we have to check distance from every corner of the closest cell, this function is faw more
332 * expensive and if Vector is known to be inside LinkedCell's domain, then SetIndexToVector() should be used.
333 * \param *x Vector with coordinates
334 * \return minimum squared distance of cell to given vector (if inside of domain, distance is 0)
335 */
336double LinkedCell::SetClosestIndexToOutsideVector(const Vector * const x) const
337{
338 for (int i=0;i<NDIM;i++) {
339 n[i] = (int)floor((x->at(i) - min[i])/RADIUS);
340 if (n[i] < 0)
341 n[i] = 0;
342 if (n[i] >= N[i])
343 n[i] = N[i]-1;
344 }
345
346 // calculate distance of cell to vector
347 double distanceSquared = 0.;
348 bool outside = true; // flag whether x is found in- or outside of LinkedCell's domain/closest cell
349 Vector corner; // current corner of closest cell
350 Vector tester; // Vector pointing from corner to center of closest cell
351 Vector Distance; // Vector from corner of closest cell to x
352
353 Vector center; // center of the closest cell
354 for (int i=0;i<NDIM;i++)
355 center[i] = min[i]+((double)n[i]+.5)*RADIUS;
356
357 int c[NDIM];
358 for (c[0]=0;c[0]<=1;c[0]++)
359 for (c[1]=0; c[1]<=1;c[1]++)
360 for (c[2]=0; c[2]<=1;c[2]++) {
361 // set up corner
362 for (int i=0;i<NDIM;i++)
363 corner[i] = min[i]+RADIUS*((double)n[i]+c[i]);
364 // set up distance vector
365 Distance = (*x) - corner;
366 const double dist = Distance.NormSquared();
367 // check whether distance is smaller
368 if (dist< distanceSquared)
369 distanceSquared = dist;
370 // check whether distance vector goes inside or outside
371 tester = center -corner;
372 if (tester.ScalarProduct(Distance) < 0)
373 outside = false;
374 }
375 return (outside ? distanceSquared : 0.);
376};
377
378/** Returns a list of all TesselPoint with distance less than \a radius to \a *Center.
379 * \param radius radius of sphere
380 * \param *center center of sphere
381 * \return list of all points inside sphere
382 */
383LinkedCell::LinkedNodes* LinkedCell::GetPointsInsideSphere(const double radius, const Vector * const center) const
384{
385 const double radiusSquared = radius*radius;
386 TesselPoint *Walker = NULL;
387 LinkedNodes *TesselList = new LinkedNodes;
388 LinkedNodes *NeighbourList = NULL;
389
390 // set index of LC to center of sphere
391 const double dist = SetClosestIndexToOutsideVector(center);
392 if (dist > 2.*radius) {
393 DoeLog(1) && (eLog()<< Verbose(1) << "Vector " << *center << " is too far away from any atom in LinkedCell's bounding box." << endl);
394 return TesselList;
395 } else
396 DoLog(1) && (Log() << Verbose(1) << "Distance of closest cell to center of sphere with radius " << radius << " is " << dist << "." << endl);
397
398 // gather all neighbours first, then look who fulfills distance criteria
399 NeighbourList = GetallNeighbours(2.*radius-dist);
400 //Log() << Verbose(1) << "I found " << NeighbourList->size() << " neighbours to check." << endl;
401 if (NeighbourList != NULL) {
402 for (LinkedNodes::const_iterator Runner = NeighbourList->begin(); Runner != NeighbourList->end(); Runner++) {
403 Walker = *Runner;
404 //Log() << Verbose(1) << "Current neighbour is at " << *Walker->node << "." << endl;
405 if ((center->DistanceSquared(*Walker->node) - radiusSquared) < MYEPSILON) {
406 TesselList->push_back(Walker);
407 }
408 }
409 delete(NeighbourList);
410 } else
411 DoeLog(2) && (eLog()<< Verbose(2) << "Around vector " << *center << " there are no atoms." << endl);
412 return TesselList;
413};
Note: See TracBrowser for help on using the repository browser.